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Figure  6-9  Lift  Coefficient  for  AR  =  2  Pitching-Plunging  Plate;  Time  History  and  vs.  Angle  of  6-10 

Attack 

Figure  6-10  Drag  Coefficient  for  AR  =  2  Pitching-Plunging  Plate;  Time  History  and  vs.  Angle  of  6-11 

Attack 

Figure  6-11  Lift  Coefficient  for  AR  =  2  Plate  in  Pure-Plunge;  Time  History  and  vs.  Angle  of  6-12 

Attack 

Figure  6-12  Drag  Coefficient  for  AR  =  2  Plate  in  Pure-Plunge;  Time  History  and  vs.  Angle  of  6-12 

Attack 
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Figure  6-13  DLR  PIV;  k  =  0.25  h  =  0.206  Pure-Plunge  of  AR  =  2  Plate,  with  a(t)  ~  Pitch-Plunge  6-13 
of2  =  0.6 

Figure  6-14  DLR  PIV;  k  =  0.125  h  =  0.412  Pure-Plunge  of  AR  =  2  Plate,  with  a(t)  ~  Pitch-Plunge  6-14 
of/l  =  0.6 

Figure  7-1  AFRL  PIV;  Re  =  40  K;  Vorticity  Contours,  k  =  3.93  Plunge;  Phases  Phi  =  0,  90,  180  7-2 

and  270;  Phase-Averaged  Measurements  and  Representative  Instantaneous 
Measurements 

Figure  7-2  AFRL  PIV;  Re  =  40  K;  Velocity  Contours,  k  =  3.93  Plunge;  Phases  Phi  =  0,  90,  180  7-2 

and  270 

Figure  7-3  AFRL  LES  Computations;  Re  =  40  K;  Vorticity  Contours,  k=3.93  Plunge;  Phases  7-3 

Phi  =  0,  90,  180  and  270;  Phase-Averaged  Measurements  and  Representative 
Instantaneous  Measurements 

Figure  7-4  AFRL  LES  Computations;  Re  =  40  K,  Velocity  Contours,  k=3.93  Plunge;  Phases  7-4 

Phi  =  0,  90,  180  and  270 

Figure  7-5  AFRL  LES  Computations;  Re  =  10  K;  Vorticity  Contours,  k  =  3.93  Plunge;  Phase  7-4 

Phi  =  0;  Phase- Averaged  and  Instantaneous;  and  Phase- Averaged  AFRL  PIV 

Figure  7-6  UM  RANS  Computations;  Re  =  40  K;  Velocity  Contours,  k  =  3.93  Plunge  7-5 

Figure  7-7  UM  RANS  Computations;  Re  =  40  K;  Vorticity  Contours,  ^  =  3.93  Plunge  7-5 

Figure  7-8  Vortex  Particle  Calculation  for  ^  =  3.93  Plunge  7-6 

Figure  7-9  Lift  Coefficient  Time  History  and  vs.  Angle  of  Attack  (right);  k  =  3.93  h  =  0.05  7-7 

Pure-Plunge  of  SD7003  Airfoil  at  Re  =  40  K  and  10  K 

Figure  7-10  Drag  Coefficient  Time  History;  k  =  3.93  h  =  0.05  Pure-Plunge  of  SD7003  Airfoil  7-7 

at  Re  =  40  K  and  10  K 

Figure  8-1  Commanded  Wing  Kinematics  Plotted  vs.  Time,  for  First  Quarter- Stroke  8-2 

Accelerating  Over  0.25  c;  Angular  Velocity  and  Acceleration  as  a  Function  of  Time 

Figure  8-2  Commanded  Wing  Kinematics  Plotted  vs.  Angular  Displacement,  for  the  First  8-2 

Quarter-Stroke;  Angular  Velocity  for  a  Wing  Accelerating  Over  the  First  0.25  c 
of  Travel 

Figure  8-3  Lift  Coefficients  for  a  Waving  Wing  Accelerating  Over  0.10  c,  0.25  c,  0.60  c  and  8-3 

All  Exponential  Velocity  Profiles 

Figure  8-4  Angular  Acceleration  for  the  Three  Velocity  Profiles:  Linear,  Sinusoidal  and  8-4 

Exponential  Over  0.10  c,  0.25  c  and  0.60  c 

Figure  8-5  Comparison  of  Unsteady  Lift  Coefficients  Observed  for  a  Waving  Wing  and  8-5 

Translating  Wing  Accelerating  Over  0.60  c 

Figure  8-6  Lift  Coefficient  as  a  Function  of  Angle  of  Attack  8-6 

Figure  8-7  Flowfield  Observed  at  Three  Points  During  the  Wing  Stroke  8-7 

Figure  8-8  Normalized  Circulation  Computed  from  PIV  Data  at  Half-Span;  AR  =  4,  a  =  25  deg  8-8 

Figure  8-9  Chordwise  View  of  Vorticity  (co/U'c)  and  A-Contours;  AR  =  4,  a  =  25  deg  8-9 

Figure  8-10  Spanwise  View  of  Vorticity  (co/U'c)  and  A-Contours;  AR  =  4,  a  =  25  deg  8-10 

Figure  8-11  Proposed  Model  of  Three  Phases  of  Flowfield  Development  for  the  Whirling  Plate  8-11 
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Unsteady  Aerodynamics  for 
Micro  Air  Vehicles 
(RTO-TR-AVT-149) 


Executive  Summary 


The  flowfield  environments  encountered  by  Micro  Air  Vehicles  (MAVs)  are  fundamentally  unsteady  - 
whether  for  fixed-wing,  rotary-wing  or  flapping-wing  configurations.  The  AVT-149  Task  Group  addressed 
fundamental  questions  in  unsteady  low  Reynolds  number  aerodynamics,  studying  computationally, 
experimentally  and  analytically  a  progression  of  abstracted  configurations,  from  2D  airfoils  in  pure-plunge 
through  3D  wing  planforms  in  plunge  and  pitch-plunge.  The  result  was  a  better  understanding  of  the 
interplay  between  the  motion,  the  flowfield,  and  the  aerodynamic  loads  history. 

Research  objectives  included: 

1)  Extension  of  prior  work  on  airfoil  laminar  boundary  layers,  separation  bubbles  and  transition  to 
unsteady  cases,  where  large  excursions  in  effective  angle  of  attack  and  production  of  concentrated 
vorticity  interacts  with  the  transition  process  to  affect  integrated  forces  and  moments.  The  result  is 
an  assessment  of  low  Reynolds  number  effects  for  MAVs,  and  survey  of  how  off-the-shelf  and 
research-type  computational  tools  can  cope  with  the  subject  flowfields. 

2)  Assessment  of  the  relation  between  vortex  shedding  (for  example,  from  airfoil  leading  edges)  and 
the  time-dependent  integrated  aerodynamic  forces  and  moments.  How  quasi-steady  is  the  force 
response,  and  how  does  it  relate  to  classical  dynamic  stall?  How  important  is  resolution  of  the 
flowfield  in  capturing  accurately  the  time  history  of  lift,  drag/thrust,  and  pitching  moment? 

3)  Assessment  of  configuration  effects;  we  compare  a  2D  airfoil  (the  SD7003),  a  2D  thin  flat  plate 
with  round  edges,  and  an  aspect  ratio  =  2  thin  flat  plate  with  round  edges. 

4)  Benchmarking  of  recently-built  facilities  and  rigs;  validation  of  recently-developed  computational 
fluid  dynamics  codes;  and  progress  in  flowfield  diagnostics  methods. 

By  comparing  results  from  14  research  groups,  the  following  conclusions  were  reached,  in  brief: 

•  Airfoils  at  moderate  Reynolds  number  (10,000  -  60,000)  and  near-stall  peak  effective  angle  of 
attack  are  strongly  sensitive  to  boundary  layer  transition  effects. 

•  On  the  other  hand,  high  effective  angle  of  attack,  well  beyond  stall,  means  that  effects  of  transition 
are  of  secondary  importance,  at  least  for  lift  production. 

•  Flat  plates  with  round  edges  have  much  lower  sensitivity  to  either  Reynolds  number  or  boundary 
layer  transition  effects. 

•  Moderate  aspect  ratio  plates  (AR  =  2)  evince  flowfield  features  irreducible  to  infinite  aspect  ratio, 
but  the  lift  and  drag  history  is  more  quasi-steady  and  more  similar  to  simple  theoretical  prediction 
than  what  one  finds  for  a  wall-to-wall  or  2D  plate. 

•  Lift  coefficient  history  is  more  quasi-steady  than  the  large  flow  separation  would  imply,  whence 
low-level  methods  still  have  good  potential  for  MAV  aerodynamics  prediction. 

•  High-frequency  low-amplitude  oscillations  are  dominated  by  non  circulatory  effects,  whence 
classical  planar- wake  models  have  good  predictive  utility. 

•  Comparison  of  translational  vs.  rotational  impulsive-start  problems  shows  no  appreciable 
difference  in  lift  history,  and  neither  case  evinces  a  stable,  attached  leading  edge  vortex. 
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Aerodynamique  instable  pour 
micro-vehicules  aeriens 
(RTO-TR-AVT-149) 


Synthese 

Les  ecoulements  rencontres  par  les  micro-vehicules  aeriens  (MAV)  sont  fondamentalement  instables  -  que 
ce  soit  pour  les  voilures  fixes,  les  voilures  toumantes  ou  les  ailes  battantes.  Le  Groupe  de  Travail  AVT-149  a 
traite  des  questions  fondamentales  sur  T  aerodynamique  instable  a  faible  nombre  de  Reynolds,  en  etudiant 
par  le  calcul,  T experimentation  et  Tanalyse  un  ensemble  de  configurations  abstraites,  allant  du  profit  2D  en 
immersion  simple  a  la  forme  en  plan  d’aile  3D  en  immersion  et  en  immersion  laterale.  II  en  est  resulte  une 
meilleure  comprehension  de  Tinteraction  entre  le  deplacement,  Tecoulement  et  les  charges  aerodynamiques. 

Les  objectifs  de  la  recherche  ont  inclus  : 

1)  Une  extension  des  travaux  prioritaires  sur  les  couches  limites  laminaires  du  profit,  les  bulles  de 
separation  et  la  transition  vers  des  cas  d’instabilite  avec  d’importantes  digressions  sur  Tangle 
d’attaque  reel  et  la  production  de  tourbillons  concentres  qui  agissent  sur  le  processus  de  transition 
affectant  les  forces  integrees  et  les  moments.  II  en  a  resulte  une  evaluation  des  effets  du  faible 
nombre  de  Reynolds  pour  les  MAV,  et  une  etude  expliquant  pourquoi  des  outils  informatiques  de 
recherche  classiques  sur  etagere  peuvent  faire  face  au  probleme  de  Tecoulement. 

2)  Une  evaluation  de  la  relation  entre  la  chute  du  vortex  (par  exemple,  a  partir  du  bord  d’attaque  du 
profit)  et  les  forces  et  les  moments  aerodynamiques  temporaires  integres.  Comment  la  reponse  de 
la  force  est  quasi  stable  et  comment  elle  correspond  a  un  decrochage  dynamique  classique  ? 
L’ importance  de  la  resolution  de  Tecoulement  dans  la  capture  precise  de  Thistorique  de  portance, 
de  trainee/poussee  et  le  moment  de  tangage  ? 

3)  Une  evaluation  des  effets  lies  a  la  configuration  ;  comparaison  du  profit  2D  (le  SD7003),  du  plan 
horizontal  fin  2D  avec  des  bords  arrondis  et  un  plan  horizontal  fin  avec  des  bords  arrondis  d’un 
allongement  geometrique  =  2. 

4)  L’etalonnage  concurrentiel  des  installations  et  des  equipements  ;  la  validation  de  codes  dynamiques 
des  fluides  informatises  recemment  developpes  ;  et  les  progres  dans  les  methodes  de  diagnostiques 
des  ecoulements. 

En  comparant  les  resultats  de  14  groupes  de  recherche,  les  conclusions  suivantes  ont  etc  obtenues, 
en  resume  : 

•  Les  profils  avec  un  nombre  de  Reynolds  modere  (10,000  -  60,000)  et  un  angle  d’attaque  reel  pres 
du  pic  de  decrochage  sont  fortement  sensibles  aux  effets  de  la  couche  limite  de  transition. 

•  D’ autre  part,  un  angle  d’attaque,  au-dela  du  decrochage,  signifie  que  les  effets  de  transition  sont 
secondaires  au  moins  en  ce  qui  concerne  la  portance. 

•  Les  plans  horizontaux  avec  des  bords  arrondis  sont  moins  sensibles  au  nombre  de  Reynolds  ou 
aux  effets  de  la  couche  limite  de  transition. 

•  Des  plans  avec  un  allongement  geometrique  modere  (AG  =  2)  presentent  des  caracteristiques 
d’ecoulement  irreductibles  a  Tallongement  geometrique  infini,  mais  Thistorique  de  la  portance  et 
de  la  trainee  est  quasi  stable  et  semblable  a  celui  de  la  prediction  theorique  simple  que  Ton  trouve 
pour  un  «  mur  contre  mur  »  ou  un  plan  2D. 
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•  L’historique  du  coefficient  de  portance  est  plus  stable  que  ce  que  la  separation  d’ecoulement  large 
devrait  impliquer,  cependant  les  methodes  a  niveau  faible  donnent  toujours  un  bon  potentiel  de 
prediction  aerodynamique  des  MAV. 

•  Les  oscillations  hautes  frequences  a  faible  amplitude  sont  dominees  par  des  effets  non  circulaires, 
cependant  les  modeles  classiques  de  sillage  planaire  sont  utiles  pour  une  bonne  prediction. 

•  La  comparaison  des  problemes  de  translation  avec  les  problemes  de  rotation  a  demarrage  impulsif 
montre  peu  de  difference  avec  I’historique  de  la  portance  et  aucun  cas  ne  montre  un  vortex  de 
bord  d’attaque  attache  et  stable. 
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Chapter  1  -  INTRODUCTION 

1.1  PHILOSOPHICAL  PRELIMINARIES 

Now  well  into  the  second  century  of  powered  flight,  we  have  much  cause  to  wonder  what  great  advances 
shall  be  bequeathed  by  our  generation,  and  whether  what  passes  for  research  in  aeronautics  is  not  an 
excess  of  the  first  syllable  and  dearth  of  the  second.  Fortunately,  complicated  problems  are  never  truly 
solved.  Perhaps  the  great  remaining  challenges  in  aeronautics,  which  are  least  solved,  lie  at  the  extremes: 
the  very  fast,  and  the  very  small  and  slow.  The  former,  in  large  measure,  are  problems  of  propulsion, 
heat  transfer  and  the  like.  But  the  latter  are  fundamental  questions  of  vorticity  transport,  and  of  the 
behaviour  of  separated  flows;  in  short,  they  are  questions  of  classical  aerodynamics,  still  fresh  even  in  our 
time,  and  still  complex. 

The  so-called  “Micro  Air  Vehicles”,  or  MAVs,  have  been  a  great  boon  to  aerospace  engineers  eager  for 
the  elusive  combination  of  relevance  to  modem  applications  and  substantive  fundamental  problems. 
MAVs  are  therefore  not  only  a  topical  application,  but  an  intriguing  problem  of  fundamental  interest  to  the 
fluid  mechanicist,  the  aeronautical  designer  and  the  biologist,  as  the  relation  between  flight  of  natural 
creatures  and  MAVs  is  not  merely  a  metaphorical  motivation,  but  is  quite  literally  tme.  The  definition  of 
MAVs  is  somewhat  amorphous,  depending  on  the  biases  of  the  definer.  Loosely  following  Pines  and 
Bohorquez  [1],  Shyy  et  al.  [2],  and  Mueller  [3],  we  can  “define”  MAVs  for  present  purposes  as  flyers  in 
the  Reynolds  number  range  of  lO"^  to  10^  based  on  the  relevant  length  scale  (typically  wing  chord)  and 
velocity  scale  (typically  flight  speed),  which  translates  into  vehicle  maximal  dimensions  on  the  order  of 
10  -  30  cm. 

MAVs  and  small  birds  share  a  direct  connection  in  size,  speed,  flight  regime  and  even  mission.  They  share 
the  same  problems  of  manoeuvring  to  avoid  obstacles  or  to  reach  points  of  interest  precisely,  to  negotiate 
gusts  and  other  ambient  flow  disturbances,  and  to  combine  lift-propulsion-control  in  a  robust,  efficient 
manner.  While  doing  so  depends  on  favourable  confluence  of  so  many  subjects,  from  flight  mechanics  to 
materials  and  actuation  and  efficient  energy  storage,  the  core  problem  for  aeronautical  conceptual  design 
and  flight  requirements  definition  is  the  relationship  between  the  aerodynamic  loads  and  the  kinematics  of 
motion  or  the  vehicle  “state”.  Motion  kinematics  may  come  from  flapping  wings  or  any  other  relative  motion 
between  the  fluid  and  the  body.  The  over-arching  condition  is  that  motion  kinematics  are  aggressive  -  fast, 
and  of  high  amplitude  -  relative  to  the  usual  aeronautical  engineering  time/length  scales  developed  for  fixed- 
wing  aircraft  in  cruise  or  gentle  manoeuvre.  Exploring  the  consequences  of  such  motion  kinematics,  in  an 
abstracted  sense,  is  the  subject  of  the  present  report. 

Unsteady  aerodynamics  in  two  and  three  dimensions  is  a  critical  area  for  scientific  understanding  of  Micro 
Air  Vehicle  flight,  for  without  it,  we  return  to  haphazard  beating  of  the  air  much  akin  to  the  first  attempts 
in  powered  flight.  Applications  of  unsteady  aerodynamics  include  perching,  gust  response,  manoeuvring 
flight,  flapping  wings,  aeroelastic  interactions  and  in  general  all  sorts  of  departures  from  the  steady  regime 
-  which,  for  small  vehicles,  includes  the  majority  of  the  flight  envelope.  For  engineering  applications, 
the  central  question  is  whether  the  aerodynamic  loads  are  essentially  quasi-steady  with  the  motion 
kinematics.  For  example  if  the  lift  and  pitch  coefficients  are  quasi-steady  with  effective  angle  of  attack, 
then  the  distinction  from  steady  aerodynamics  is  trivial.  The  secondary  question  is  the  extent  to  which  the 
relative  motion  between  body  and  fluid  forces  the  aerodynamic  response  such  that  the  usual  variables  in 
aerodynamics  -  Reynolds  number,  boundary  layer  transition,  roughness  and  leading  edge  curvature, 
and  geometry  in  general  -  become  relatively  unimportant;  or,  if  these  various  factors  require  consideration 
just  as  they  do  in  steady  aerodynamic  problems. 

We  begin  by  reviewing  the  classical  problem  of  pitching  and  plunging  airfoils,  and  then  consider  its  relation 
to  the  aerodynamics  of  flight  in  nature.  Next,  we  propose  a  set  of  canonical  problems  for  experiments  and 
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computations  in  low  Reynolds  number  unsteady  aerodynamics.  After  describing  the  various  facilities  and 
computational  methods  used  to  obtain  the  results  for  the  canonical  problems,  the  bulk  of  this  report  recites 
the  various  results  case-by-case,  drawing  hopefully  insightful  inferences  for  each.  We  then  consider  related 
problems  not  covered  in  the  main  under  the  canonical  case,  and  conclude  with  recommendations  for 
generalizing  the  canonical  problems  and  for  other  points  of  departure. 


1.2  RESUME  OF  PRIOR  RESEARCH  IN  2D  PITCHING  AND  PLUNGING 

In  two  dimensions  the  possible  motions  are  rotation,  or  pitch;  up  and  down  translation,  or  plunge  - 
sometimes  called  heave;  and  fore  and  aft  translation,  or  surge.  The  latter  we  largely  ignore  in  the  present 
work,  because  of  the  complexity  of  realizing  this  degree  of  freedom  in  experiments,  and  because  it  is 
usually  absent  in  classical  pitch-plunge  problems.  Some  of  the  experimental  setups  used  in  the  present 
study  can,  however,  accommodate  all  3  longitudinal  degrees  of  freedom. 

For  all  of  the  geometric  and  kinematic  complexity  of  the  MAV  problem,  especially  for  flapping  wings, 
2D  problems  remain  important.  Simply  put,  one  has  to  begin  somewhere,  especially  if  the  context  is 
passage  from  the  steady  case.  Classical  airfoil  problems,  steady  and  unsteady,  are  framed  in  2D,  and  the 
present  approach  is  to  begin  with  2D  aerodynamics  in  the  longitudinal  plane  (lift,  pitch,  drag/thrust), 
then  to  extend  to  longitudinal-type  motions  of  a  3D  planform,  and  finally  to  consider  non-longitudinal 
motions  such  as  rotations. 

Perhaps  the  first  appearance  of  pitching/plunging  airfoil  physics  was  in  applications  of  fixed-wing  aircraft 
aeroelasticity.  Prior  to  current  interest  in  flapping  wing  aerodynamics,  dynamic  stall  of  helicopter  blades  was 
the  main  application  for  unsteady  airfoils  that  included  significant  motion  amplitude  as  well  as  frequency, 
resulting  in  the  so-called  dynamic  stall  [5].  Small  reduced  frequencies,  meanwhile,  traditionally  fell  under 
the  more  usual  problem  of  airplane  flight  dynamics,  in  the  sense  of  phugoid  and  short-period-mode 
oscillations  [4].  But  even  the  short-period  mode  is  low  frequency  in  the  sense  of  unsteady  aerodynamics, 
and  almost  without  question  is  amenable  to  quasi-steady  treatment.  Helicopter  blade  stall  problems  are 
another  matter,  introducing  at  least  two  fundamental  problems.  One  is  how  to  account  for  shedding  of 
vorticity  into  the  wake  whenever  there  is  a  change  in  effective  angle  of  attack,  either  due  to  geometric  pitch 
change  or  to  plunge.  This  is  an  essentially  inviscid  problem.  The  second  problem  is  how  to  account  for  flow 
separation,  of  which  the  most  celebrated  example  is  the  leading  edge  vortex  (LEV),  whose  formation, 
residency  time,  shedding  and  downstream  convection  form  the  core  problem  of  dynamic  stall.  McCroskey 
et  al.  [5]  pointed  out  that  as  this  vortex  passes  over  the  airfoil  surface,  it  significantly  changes  the  chordwise 
pressure  distribution  and  produces  transient  forces  and  moments  fundamentally  different  from  those  in  static 
stall.  Comprehensive  reviews  of  dynamic  stall  are  given  by  McCroskey  [6],  Carr  [7],  Liiva  [8]  and  Carr  and 
McCroskey  [9].  These  works  inform  the  origin  of  MAV-related  pitch-plunge  problems,  but  now  the 
Reynolds  number  is  a  factor  of  0(100)  smaller  than  in  helicopter  applications,  and  in  the  past  20-30  years 
this  has  led  to  a  spate  of  investigations  on  what  happens  with  flow  separation  due  to  airfoil  motions  at  low 
Reynolds  number.  The  motivation  for  these  works  is  not  necessarily  MAVs,  as  the  MAV-relevant  Reynolds 
number  just  happens  to  be  the  operational  Reynolds  number  in  many  small-scale  facilities,  especially  in 
water  tunnels  and  water  tow  tanks,  and  therefore  many  fundamental  studies  not  motivated  by  MAVs 
nevertheless  speak  to  MAV-relevant  results.  A  comprehensive  review  is  beyond  the  scope  of  the  present 
work,  but  our  purpose  is  to  systematize  the  available  results  in  a  new,  clean-sheet  effort  that  attempts  to 
benchmark  the  various  experimental  and  computational  approaches,  opening  vistas  for  parameter  studies  for 
targeted  applications. 

By  far  the  most  common  configuration  in  2D  problems  is  the  venerable  NACA  0012  airfoil  -  which  is 
perhaps  regrettable,  because  the  NACA  0012  performs  poorly  at  MAV  Reynolds  numbers.  Ohmi  et  al. 
[10,11]  considered  a  sinusoidally -pitching  NACA0012  airfoil  with  impulsive-start  to  steady  translation  in 
a  water  tow  tank.  Both  the  reduced  frequency  and  amplitude  were  found  to  affect  the  structure  of  the 
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vortex  wake.  Visbal  and  Shang  [12]  studied  the  high-rate  high-angle  of  attack  linear  pitch-ramp  of  a 
NACA0015  airfoil  at  Re  =  lO"^  by  solving  the  full  2D  Navier-Stokes  equations,  finding  that  lags  between 
evolution  of  leading-edge  flow  separation  and  the  airfoil  motion  kinematics  should  increase  with 
increasing  reduced  frequency.  Ghosh  Choudhuri  and  Knight  [13]  examined  the  effects  of  compressibility, 
pitch  rate,  and  Reynolds  number  on  the  initial  stages  of  2D  unsteady  separation  of  laminar  subsonic  flow 
over  a  pitching  airfoil  in  the  Reynolds  number  ranging  from  lO"^  to  10^,  finding  that  increasing  the 
Reynolds  number  hastens  the  appearance  of  the  primary  re-circulating  region. 

The  aforementioned  studies  focus  mostly  on  transients  following  the  initiation  of  the  airfoil  motion  from 
rest.  Others  considered  the  periodic  or  phase-averaged  behaviour  of  pitch-plunge  motions  after  initial 
transients  have  relaxed,  typically  with  a  focus  on  motion  kinematics  for  optimal  thrust  efficiency.  Platzer 
and  Jones  [14]  discussed  theoretical  prediction  of  thrust  efficiency  vs.  flow  visualization  and  thrust 
measurements  for  an  airfoil  in  pure-plunge,  over  a  range  of  reduced  frequencies  and  reduced  amplitudes. 
Koochesfahani  studied  high-frequency  low-amplitude  airfoil  pitch  oscillations  with  quantitative  and 
qualitative  visualization,  identifying  vortex  shedding  patterns  vs.  motion  kinematic  parameters  [15]. 
McAlister  and  Carr  [16],  and  Walker  et  al.  [17]  examined  the  interplay  between  leading-edge  and  trailing- 
edge  vortices.  Anderson  et  al.  [18]  and  Triantafyllou  et  al.  [19]  used  particle  image  velocimetry  and  direct 
force  measurement  to  study  the  combined  pitch-plunge  parameter  space  for  propulsive  efficiency 
optimization,  conducting  a  large  parameter  study,  albeit  force  measurements  and  flowfield  data  were 
obtained  at  quite  disparate  Reynolds  numbers.  Young  and  Lai  [20]  used  a  2D  Reynolds-Averaged  Navier- 
Stokes  (RANS)  approach  to  study  the  frequency-amplitude  parameter  space  for  optimal  thrust  efficiency 
for  a  plunging  NACA  0012  airfoil  at  Re  =  0(10"^),  and  Jones  et  al.  [21]  demonstrated  good  agreement  in 
wake  vortex  structure  between  dye  injection  at  the  airfoil  trailing  edge  in  a  water  tunnel,  2D  laminar 
Navier-Stokes  computation,  and  a  2D  vortex-particle  method.  Lian  and  Shyy  [22]  showed  similar 
agreement  using  a  RANS  solver  run  fully-turbulent.  Here  we  seek  to  extend  these  results  to  more  complex 
kinematics  of  motion,  to  a  range  of  Reynolds  numbers  bracketing  transition-dominated  effects  evinced  in 
the  static  measurements  and  in  lower-amplitude  pitch-plunge  oscillations  [23],  and  to  include  effects  of 
leading-edge  separation. 

One  important  issue  in  periodic  oscillatory  airfoil  flows  is  the  lag  between  the  aerodynamic  response  and 
the  airfoil  motion  kinematics  in  2-degree-of-freedom  pitch-plunge.  As  a  natural  extension  of  the  strictly 
quasi-steady  model,  one  seeks  to  construct  an  explicit  relation  of  the  lag  of  putatively  sinusoidal  force 
response  to  sinusoidal  motion  kinematics,  as  a  function  of  reduced  frequency,  of  amplitudes  of  pitch  and 
plunge,  of  phase  difference  between  pitch  and  plunge,  of  Reynolds  number  and  so  forth.  This  could  then 
form  a  model  for  the  lift  response  to  more  general  motions  and  in  more  general  configurations,  serving  as 
an  engineering  tool  for  parameter  studies. 


1.3  HIERARCHY  OF  ANALYSIS  FOR  UNSTEADY  AIRFOIL  AND  WING 
AERODYNAMICS 

1.3.1  Classical  Methods  and  Closed-Form  Solutions 

Analytical  methods  in  aerodynamics  are  most  successful  for  small  disturbances,  and  in  unsteady 
aerodynamics  this  means  small  amplitudes  of  motion,  though  not  necessarily  low  frequencies.  Analytical 
methods  for  the  aerodynamics  of  moving  airfoils  and  wings  date  back  to  the  period  of  1920  -  1940  when 
increasing  aircraft  flight  speeds  pushed  the  need  to  estimate  and  to  analyse  flutter  instabilities  of  flight 
surfaces.  This  therefore  predates  the  helicopter-motivated  problems  cited  earlier.  The  early  work  of 
Bimbaum  [24]  solved  the  unsteady,  inviscid  flow  problem  of  flat  plates  for  harmonic  low-frequency  pitch 
and  plunge  motions.  Improved  ansatz  functions  along  the  frequency  were  found  independently  by  Kiissner 
[25]  and  Theodorsen  [26],  and  these  yielded  closed-form  solutions  for  the  aerodynamic  forces  of 
harmonically  oscillating  lifting  airfoils  at  small  amplitude  but  with  no  restriction  on  frequency.  While 
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starting  and  stopping  vortices  in  the  airfoil  wake  are  modelled,  the  wake  shape  is  assumed  to  be  a  straight 
line.  Theodorsen’s  solution  was  used  by  Garrick  [27]  to  compute  drag  or  thrust  as  well.  Along  similar 
lines  Wagner  [28]  derived  the  lift  function  for  the  problem  of  impulsive  start,  which  is  more  suited  for 
hover  problems. 

These  theories  may  be  used  for  rapid  calculations  for  parameter  studies,  yielding  reasonable  estimates  of 
thrust  and  propulsive  efficiency  for  2D  problems,  and  they  can  be  used  as  the  2D  element  of  strip  theories 
for  wing  flows.  But  there  are  two  key  obstacles.  First,  what  is  conventionally  lumped  as  “viscous  effects” 
must  be  small,  for  example  limited  to  skin  friction  or  laminar  separations  closed  in  the  time-averaged 
sense.  Second,  3D  effects  must  be  small,  in  that  spanwise  variations  can  be  accounted  by  a  local  effective 
angle  of  attack.  Both  assumptions  in  general  fail  for  flapping  wing  aerodynamics  at  large  effective  angles 
of  attack  or  for  small  aspect  ratios.  Thus  there  is  a  need  to  extend  the  classical  methods  to  applications 
such  as  insect  flight  and  the  flapping-wing  MAVs,  ideally  retaining  the  simplicity  of  closed-form 
solutions. 


1.3. 1.1  Further  Remarks  on  Theodorsen’s  Method 

In  Theodorsen’s  model  [26],  for  harmonic  motions  in  pitch  and  plunge,  the  lift  coefficient  is  given  by: 
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The  pitch  and  plunge  are  input  as  complex  exponentials,  a{t)  =  aQ+ and  h(t)  = 

The  phase  lead  of  pitch  vs.  plunge,  in  fractions  of  motion  period,  is  cp.  The  reduced  frequency,  k,  is  defined 
as  ^  =  cocIlU^  =  nfc/Uao.  The  first  term  in  Eqn.  1  is  the  steady-state  lift.  The  second  term  is  the  “apparent 
mass”  or  non-circulatory  lift,  due  to  acceleration  effects.  It  is  progressively  larger  for  larger  k.  The  third 
term  models  circulatory  effects.  C(k)  is  the  complex-valued  “Theodorsen  function”,  with  magnitude  <  1. 
It  is  essentially  a  transfer  function,  accounting  for  attenuation  of  lift  amplitude  and  phase  lag  in  lift  response, 
from  its  real  and  imaginary  parts,  respectively.  Setting  C(k)  =  1  (k  =  0)  and  ignoring  non-circulatory  effects 
recovers  quasi-steady  thin  airfoil  theory.  The  non-circulatory  term  follows  instantaneously  the  kinematics  of 
motion,  but  evolution  of  vorticity  in  the  wake  yields  phase  lag  relative  to  the  kinematics  of  airfoil  motion  in 
the  circulatory  term,  which  is  predicted  to  peak  for  ^  ~  0.3. 

Equation  (1)  assumes  a  planar  wake  and  a  trailing-edge  Kutta  condition,  thus  excluding  wake  rollup  and 
vortex  streets,  vortex  shedding,  convection  of  large  separations  over  the  airfoil,  open  separations,  large 
laminar  separation  bubbles,  leading  edge  and  trailing  edge  vortices  and  so  forth.  But  the  simplicity  of  such  a 
model  obviously  makes  it  attractive  as  an  engineering  tool.  Thus,  we  consider  to  what  extent  this  model  is 
useful  for  kinematics  where  the  flow  separation  in  portions  of  the  airfoil  oscillation  ranges  from  moderate  to 
severe,  and  compare  the  Theodorsen  prediction  with  that  of  a  range  of  computations  and  experiments. 


1.3. 1.2  Modern  Extensions  of  Classical  Methods 

Aerodynamic  analysis  of  insect  flight  indicates  that  the  classical  airfoil  theories  together  with  present 
knowledge  in  airfoil  stall  result  in  much  too  low  lift  coefficients  as  compared  to  those  observed  with  natural 
flyers.  It  was  concluded  that  wing- wake  interactions  and  leading  edge  separation  contribute  significantly  to 
the  overall  forces  and  semi-empirical  models  for  these  effects  were  recently  developed,  i.e.  by  Penderson  and 
Zbikowski  [29]. 

Extension  of  airfoil  theories  to  thick  airfoils  and  modelling  the  effect  of  deforming  wakes  is  possible  with 
panel  methods,  which  approximate  the  required  vortex  distributions  along  the  airfoil  and  in  the  wake  using 
discrete  surface  panels  [30,21].  However,  the  numerical  discretization  and  iterative  determination  of  the 
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wake  require  higher  computation  efforts,  and  the  robustness  of  correctly  capturing  separation,  shedding, 
reattachment  and  other  vortical  effects  remains  a  subject  of  controversy.  More  details  on  airfoil  theory 
assessment  for  predicting  flapping  wing  propulsion  are  available  from  a  recent  review  paper  of  Platzer 
[31].  In  sum,  aerodynamic  designers  of  today  lack  a  robust  method  that  would  allow  one  to  analyse  the 
dynamics  of  inviscid  and  viscous  flowfields  of  moving  airfoils  similarly  to  the  state  of  the  art  with  the 
XFOIL  code  for  steady  flows. 

1.3.2  Viscous  Solutions  and  Modelling 

Perhaps  the  effects  of  boundary  layers  and  the  development  of  complex  flow  separations  from  edges  and 
smooth  parts  of  aerodynamic  surfaces  can  only  be  computed  by  solving  the  non-linear  equations  of  motion 
for  viscous  flows.  This  is  a  rather  straight  forward  task  for  wing  chord  Reynolds  numbers  below  values 
around  104  where  the  flow  around  the  wing  is  expected  to  be  laminar,  under  normal  circumstances; 
but  even  then,  computations  are  too  large  to  support  extensive  parameter  studies,  and  a  modelling  approach 
remains  imperative.  At  larger  Reynolds  numbers  flow  instabilities  present  in  adverse  pressure  boundary 
layers  and  in  separated  flowfields  will  eventually  transition  into  turbulence,  thereby  introducing  strong 
momentum  transport  normal  to  the  mean  streamwise  flow  directions.  The  direct  numerical  simulation  of  the 
flow  instabilities  and  their  break  down  to  turbulence  as  part  of  solving  the  equation  of  motion  of  the  unsteady 
flowfield  is  very  time  consuming  since  this  task  requires  the  numerical  resolution  of  broad  ranges  of  scales  in 
space  in  time.  It  is  no  surprise  that  this  has  only  been  accomplished  for  some  limited  research  flows  [32]  and 
there  is  no  hope  of  using  this  approach  routinely  for  exploring  the  design  space  of  practical  flapping  wings 
within  the  mid-term  future. 

Hence  the  development  of  suited  engineering  models  for  transition  and  turbulence  is  of  current  interest  to 
achieve  reasonable  prediction  capabilities  at  affordable  computation  costs.  Turbulence  models  compute 
turbulent  transport  in  terms  of  Reynolds  stresses  that  appear  as  new  unknown  after  ensemble  averaging  the 
unsteady  equations  of  motion.  The  models  consist  of  one  or  more  transport  equations  that  determine  the 
physical  behaviour  of  turbulence  quantities  and  the  needed  Reynolds  stresses  are  derived  from  these. 
The  overall  approach  is  labelled  URANS.  The  conceptually  quite  simple  application  is  then  to  use  the 
chosen  model  of  turbulence  along  the  complete  wing  surface.  However,  this  can  result  in  a  simulation  of 
fully  turbulent  flow  in  regions  where  laminar  boundary  layers  are  observed  in  experiments  and  hence, 
physical  flow  separations  may  be  suppressed  which  would  otherwise  generate  significant  contributions  to 
the  aerodynamic  forces  and  moments. 

The  instabilities  of  the  laminar  boundary  layers  that  lead  to  transition  on  airfoils  and  wings  exhibit  a 
convective  character.  It  is  then  straight  forward  to  model  transition  to  turbulence  by  assuming  an 
additional  transport  equation  for  an  intermittency  parameter  that  discriminates  between  laminar  and 
turbulent  boundary  layer  flow  regions.  This  approach  is  followed  by  Menter  et  al.  [33],  but  it  has  not  yet 
been  validated  for  Reynolds  numbers  below  10^5.  Transition  on  airfoils  and  wings  at  low  Reynolds 
numbers  depends  on  the  growth  of  harmonic  modes  in  the  laminar  boundary  layer  in  many  cases. 
This  observation  motivates  the  so-called  e^N-approach  where  the  point  of  transition  onset  is  correlated 
with  computed  amplitude  ratios  of  the  primary  boundary  layer  instabilities.  This  transition  method  has 
been  successfully  coupled  to  URANS  solvers  in  steady  and  unsteady  formulations  [34,35]  and  validation 
results  for  unsteady  low-Reynolds-number  airfoil  flows  with  moderate  flow  separations  seem  promising. 

1.3.2.1  3D  Effects 

Classical  theories  for  computing  the  effect  of  finite  span  of  flapping  wings  may  be  divided  into  unsteady 
wing  theory  for  forward  flight  and  blade  element  theory  that  is  more  suited  for  hover  and  low-speed  flight. 
While  some  of  the  potential  flow  methods  for  wings  developed  for  flutter  in  forward  flight  are  only  valid 
for  small  geometric  amplitudes  others  were  specifically  developed  for  flapping  wings.  Philips  et  al.  [36] 
developed  a  lifting  line  theory  with  a  suited  model  of  axial  and  transverse  vortices  in  the  wake.  Hall  et  al. 
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[37]  presented  an  unsteady  vortex  lattice  method  for  more  general  wing  planforms.  These  approaches  are 
used  to  assess  propulsive  efficiencies  of  wings  in  forward  flight. 

Blade  element  theory,  on  the  other  hand,  can  model  both  steady  and  unsteady  flow  effects  along  the  3D 
wing  in  hover  and  many  variants  of  this  approach  appeared  over  the  last  25  years  [38,39,40,29],  using  2D 
quasi-steady  force  coefficients  or  unsteady  2D  theories  as  elementary  data  source.  The  use  of  certain  flow 
models  in  these  wing  theories  has  generally  been  validated  along  a  small  amount  of  force  data  available 
from  dedicated  experiments  for  specific  motion  kinematics.  This  appears  as  a  major  weakness  today,  as  a 
truly  predictive  capability  for  complex  flows  can  likely  not  be  obtained  by  this  approach. 


1.3.2.2  Importance  in  Capturing  Transition  for  Low  Reynolds  Number  Unsteady  Aerodynamics 
Computations 

While  it  is  well  known  that  laminar-turbulent  transition  in  boundary  layers  has  a  large  impact  on  airfoil 
drag  at  moderate  to  high  Reynolds  numbers,  but  the  impact  in  low-Reynolds-number  airfoils  is  different. 
The  importance  of  the  transition  at  low  Reynolds  numbers  stems  from  the  observation  that  laminar 
boundary  layers  separate  more  early  under  adverse  pressure  gradients  than  their  turbulent  counterparts. 
Low-Reynolds-number  airfoils  exhibit  therefore  three  important  flow  phenomena  that  depend  on  the  angle 
of  attack  and  possibly  on  its  unsteady  variation.  These  are  laminar  separation  bubbles  (LSB),  laminar 
airfoil  stall  and  dynamic  airfoil  stall. 

The  LSB  denotes  a  flow  where  a  laminar  separation  takes  place  caused  by  an  adverse  pressure  gradient 
along  the  aerodynamic  surface.  Small  harmonic  disturbances  present  in  the  laminar  boundary  layer  are 
strongly  amplified  in  the  shear  layer  of  the  separated  flow  [41]  and  rapid  transition  to  turbulence  takes 
place.  The  turbulence,  in  turn,  creates  a  large  momentum  transport  normal  to  the  shear  layer  so  that  the 
flow  reattaches  to  the  surface  and  a  closed  bubble  is  formed  in  the  time-averaged  mean.  LSBs  can  create 
additional  pressure  drag  as  they  displace  the  outer,  inviscid  flow. 

A  more  important  aerodynamic  effect  occurs  for  steady-state  onset  flow  if  the  transition  process  in  the 
separated  shear  layer  is  relatively  slow  and  the  adverse  pressure  gradient  is  strong.  Then  turbulent 
momentum  transport  is  not  sufficient  to  close  the  bubble  and  a  large  separation  occurs  that  extends  right  to 
the  trailing  edge.  This  causes  a  sudden  loss  of  lift  and  a  strong  increase  of  drag  along  with  significant 
hysteresis  effects  of  force  coefficients  with  varying  angle  of  attack.  Note  that  the  separated  flow  region  is 
characterized  by  continuous  vortex  shedding  because  of  the  unstable  behaviour  of  the  shear  layers  present. 

For  large  unsteady  increases  of  the  airfoil  angle  of  attack  beyond  the  steady  stall  limits  the  laminar  flow 
separation  at  the  leading  edge  may  develop  into  a  strong  and  well  organized  vortex.  This  vortex  induces 
significant  suction  on  the  airfoil  surface  with  an  overshoot  of  lift  and  usually  significant  variations  of  the 
pitching  moment.  Airfoil  drag  is  also  increased  by  this  dynamic  stall  effect.  As  a  consequence  it  may  be 
assumed  that  truly  predictive  computations  should  take  into  account  transition  in  laminar  boundary  layers 
and  particular  within  laminar  separations  since  steady  and  dynamic  airfoils  stall  depends  on  the  boundary 
layer  state. 


1.4  THE  DIFFERENT  KINDS  OF  FLAPPING  IN  NATURE 

A  large  but  not  exclusive  motivation  for  our  work  is  the  engineering  abstractions  of  flapping-flight  in 
nature,  both  as  routes  for  rational  application  of  bio-inspiration  in  MAV  design,  and  for  first-principles 
understanding  of  the  aerodynamics  of  natural  flyers.  Obviously  analysis  of  the  full  spectrum  of  natural 
flight  is  hubristically  ambitious  and  is  not  attempted  here.  It  is  however  worthwhile  to  attempt  to  place  the 
present  work  in  context.  Thus  we  provide  a  very  brief  and  somewhat  superficial  summary  of  the  different 
flapping  regimes  found  in  nature,  at  least  in  terms  of  their  aerodynamic  context. 
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The  method  of  flight  best  understood  to  date  is  flxed  wing  flight,  where  lift  is  obtained  by  mainly  attached 
flow  over  a  nominally  rigid  wing.  Many  birds  are  seen  to  spend  a  lot  of  their  air-time  in  a  similar  flight 
regime,  soaring  in  thermals  or  gliding  some  distance  between  bursts  of  flapping.  However,  we  can  also 
observe  that  steady  gliding  flight  becomes  increasingly  uncommon  with  reducing  size  (compare  insects 
with  large  birds),  and  the  first  question  to  ask  is  whether  there  are  fundamental  constraints  that  prevent 
gliding  flight  at  small  scales. 

Assuming  steady  flight,  lift  equals  weight,  and  we  can  use  observations  from  zoologists  to  estimate  the 
scaling  of  wing  performance  of  birds  and  insects.  As  with  airplanes,  ^  puS  Q  =  mg,  where  p  is  the  air 

density,  U  is  the  forward  (gliding)  velocity,  S  is  the  wing  area.  Cl  is  the  wing  lift  coefficient,  mg  is  the 
total  weight.  Figure  1-1  shows  data  collected  for  a  very  large  range  of  birds,  comparing  typical  wing  areas 
with  total  mass.  Here  we  see  that  for  most  birds  (except  the  smallest  hummingbirds)  there  is  a  clear 
relationship  between  wing  area  and  body  mass  which  can  be  expressed  by  a  best  fit  as  iS  oc  . 


CONSEQUENCES  OF  SIZE  DIFFERENCES 


Figure  1-1:  Scaling  of  Bird  Wing  Area  vs.  Body  Mass  -  from  McMahon  and  Bonner  [42]. 


From  the  above  we  can  then  write  for  the  gliding  speed: 


U  = 


(2) 
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If,  for  the  moment,  we  assume  that  gliding  flight  for  all  birds  is  achieved  with  the  same  wing  lift 
coefficient,  then  the  above  would  imply  (making  use  of  the  approximate  proportionality  between  wing 
area  and  that  the  relationship  between  ‘cruise’  speed  and  body  mass  is  U .  Fortunately, 

zoologists  have  observed  bird  and  insect  flight  in  nature.  Figure  1-2  shows  typical  measured  forward 
velocities  for  a  large  number  of  animals  and  aircraft.  Again  it  can  be  seen  that  the  data  follows  a  clear 
trend  (albeit  with  fairly  large  departures  for  individual  examples).  The  best  fit  line  is  U which  is 
surprisingly  close  to  the  expected  result  from  the  above. 
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Figure  1-2:  Weight  and  Fiight  Speed  for  Animais  and  Aircraft  -  from  Tennekes  [45]. 
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From  this  simple  discussion  we  can  conclude  that  for  gliding  flight  in  nature  (and  aeronautics)  the  wing  lift 
coefficient  is  almost  completely  independent  of  vehicle/animal  size.  The  wing  on  a  small  bird  must  generate 
the  same  overall  lift  coefficient  as  that  on  a  large  aircraft.  In  fact,  closer  examination  would  suggest  that  lift 
coefficients  might  even  slightly  increase  as  size  reduces  (although  at  small  sizes  gliding  flight  is  increasingly 
rare  making  this  crude  analysis  invalid).  Others,  such  as  Ellington  [38],  have  analysed  hovering  flight  and 
concluded  the  required  lift  coefficients  are  at  least  as  great  as  these  observed  in  gliding  flight  at  larger  scales. 

Therefore,  the  wing  sections  employed  in  any  aircraft  need  to  be  capable  of  achieving  similar  lift  coefficients 
regardless  of  scale,  and  thus  Reynolds  number.  This,  however,  is  not  at  all  straightforward. 

Figure  1-3  shows  the  result  of  a  literature  survey  of  airfoil  performance  across  a  wide  range  of  Reynolds 
numbers.  Here,  it  can  be  seen  that  there  is  a  critical  range  between  roughly  Re  =  10000  and  Re  =  100000, 
where  maximum  airfoil  lift  coefficients  experience  a  clear  drop  as  a  result  of  laminar  separation. 
This  critical  range  is  roughly  in  the  regime  where  small  to  medium  sized  birds  operate.  The  performance 
of  wings  in  this  critical  range  can  be  improved  by  employing  additional  transition  mechanisms,  such  as 
roughness  or  leading  edge  devices,  and  such  techniques  can  indeed  be  found  in  nature  [43,44].  However, 
while  the  useful  range  can  be  extended  down  towards  smaller  Reynolds  numbers  with  such  flow  control, 
this  does  not  remove  the  fundamental  problem  that  below  the  critical  Re  regime  steady  airfoil  performance 
is  insufficient  for  gliding  flight. 


There  are,  however,  some  data  included  in  Figure  1-3  which  show  significantly  higher  lift  coefficients  in  the 
sub-critical  Re  regime.  These  are  found  in  experiments  that  either  feature  unsteadiness  (stopping/starting 
wings)  or  three-dimensionality  (rotating  wings).  This  suggests  that  for  smaller  sized  vehicles,  unsteady  lift 
mechanism  (or  at  least  three-dimensional  ones)  are  necessary. 
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Depending  on  the  size  of  the  animal  and  the  type  of  flight  we  can  now  observe  several  types  of  flapping 
wing  flight  in  nature  which  will  be  discussed  in  turn  below. 

1.4.1  Large  Birds  in  Forward  Flight 

For  flight  in  this  regime  the  wing  size  and  flight  velocity  are  large  enough  for  airfoil  Reynolds  numbers  to 
be  in  excess  of  the  critical  range.  Here,  flapping  is  only  necessary  to  generate  thrust,  and  since  many  birds 
are  highly  efficient  (in  terms  of  lift/drag)  this  is  not  a  large  force. 

To  generate  thrust  birds  make  use  of  the  fact  that  a  vertical  (plunging)  motion  of  the  wing  -  a  vertical 
stroke  plane  -  rotates  the  effective  angle  of  attack  vector  as  seen  in  Figure  1-4.  During  a  downstroke  this 
generates  thrust.  In  contrast,  the  upstroke  is  relatively  passive  -  the  wing  is  rotated  to  a  near  zero  effective 
angle  of  attack. 


Typically  (because  of  the  low  thrust  requirement)  only  the  outer  wing  is  used  for  thrust  generation  and  the 
inner  wing  flaps  much  less.  This  is  shown  in  Figure  1-5.  The  relatively  small  plunge  experienced  by  the 
inner  wing  ensures  that  it  continues  to  generate  lift  through  the  cycle  and  there  is  only  a  minor  variation  of 
the  lift  vector.  By  contrast  the  outer  wing,  which  undergoes  a  greater  amplitude  plunging  motion  (together 
with  pitch  variation)  produces  the  thrust  and  lift  during  the  downstroke  only.  Flapping  frequencies  are 
typically  quite  low  and  the  flowfield  is  effectively  quasi-steady  (and  attached)  throughout.  However, 
the  forces  experienced  by  the  wing  are  not  quasi-steady  as  unsteady  effects  already  have  an  impact, 
for  example  through  the  time- varying  vorticity  in  the  wake. 


Figure  1-5:  Comparison  of  inner  and  Outer  Wing  Motion  for  Condor  in  Fast  Fiapping  Fiight  [46]. 
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1.4.2  Medium  Birds  in  Forward  Flight,  Hovering  Flight  (Large  and  Medium) 

At  lower  flight  speeds,  for  example  before  landing  or  after  take-off,  there  is  insufficient  forward  velocity 
to  produce  the  required  lift.  Medium  sized  birds  are  also  approaching  the  critical  Reynolds  regime  making 
lift  generation  more  difficult.  In  such  a  situations  wing  motion  is  employed  to  achieve  greater  effective 
velocity.  The  extreme  case  is  during  hover  when  all  the  effective  free-stream  velocity  is  generated  by  the 
wing  motion.  For  flight  in  this  category  the  downstroke  is  generally  sufficient  to  produce  the  necessary 
forces  -  consistent  with  the  observation  that  downstroke-enabling  vs.  upstroke-enabling  flight  muscles  are 
highly  unequal  in  strength. 

Figure  1-6  shows  a  pigeon  and  a  bat  in  hover.  Here  the  stroke  plane  is  rotated  to  produce  a  net  force  acting 
in  vertical  direction.  During  the  upstroke  the  wings  are  folded  for  minimum  affective  area  -  the  upstroke  is 
effectively  passive.  Note  that  both  lift  and  drag  (relative  to  the  effective  local  flow  angle)  contribute  to  the 
force  production. 


Figure  1-6:  Hovering  Fiight  in  Pigeons  and  Bats  -  after  Brown  [47]. 
Lift  is  generated  during  the  downstroke  oniy,  the  upstroke  is  passive. 
The  stroke  plane  is  rotated  to  give  an  effective  vertical  force  vector. 


In  this  regime,  the  combination  of  wing  size  and  speed  of  motion  is  sufficient  to  ensure  that  the  flow  remains 
above  the  critical  Reynolds  number  range.  This  means  that  lift  is  primarily  generated  in  the  traditional 
manner  (that  is,  attached  and  quasi-steady  airfoil  lift)  during  the  down-stroke.  However,  extreme  manoeuvres 
are  likely  to  require  genuine  unsteady  and  separated  flow  aerodynamics. 

1.4.3  Small  Birds  and  Insects  in  Forward  Flight  or  Hover 

When  the  wing  size  reduces  even  further  the  flight  speed  also  reduces.  To  some  extent  this  is  offset  by  an 
increase  in  flapping  frequency,  however,  this  is  no  longer  sufficient  to  avoid  the  drop  in  lift  experienced 
due  to  two  factors:  the  Reynolds  number  is  now  in  or  below  the  critical  range;  and  the  maximum  effective 
velocity  relative  to  the  wing  is  reduced.  In  this  regime  the  downstroke  alone  is  insufficient  to  generate  the 
necessary  lift  (and  wing  muscles  are  more  equal  in  strength).  As  seen  in  Figure  1-7,  showing  a  hovering 
hummingbird,  the  stroke  plane  is  now  almost  horizontal  and  the  wing  is  pitched  at  either  end  of  the  stroke 
to  give  an  effective  positive  angle  of  attack  during  both  the  up-  and  downstroke.  Effectively  the  wing  now 
works  in  a  propeller-type  motion,  similar  to  a  helicopter.  As  a  result  of  the  horizontal  stroke -plane, 
the  total  force  vector  is  no  longer  vertical  during  each  up/downstroke.  This  is  because  the  drag  component 
introduces  a  horizontal  force.  Because  of  the  symmetry  of  the  up/downstroke  motion,  the  cycle-averaged 
horizontal  force  cancels  out. 
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Figure  1-7:  Hummingbird  in  Hover  -  after  Stoipe  and  Zimmer  [48]. 

Similar  wing  kinematics  are  also  applied  by  insects.  For  forward  flight  the  stroke  plane  is  rotated  from  the 
horizontal  to  achieve  a  thrust  component.  Gliding  flight  is  virtually  non-existent  in  this  regime  as  a  result 
of  the  poor  lift  coefficients  available  at  low  Reynolds  numbers. 

Even  at  high  flapping  frequencies,  the  effective  Reynolds  number  is  too  low  to  ensure  good  lift  coefficients 
and  in  this  regime  unsteady  and  three-dimensional  effects,  such  as  the  leading  edge  vortex,  play  a  major  role 
in  lift  generation. 

In  surveying  the  various  regimes  of  flapping-wing  flight  from  an  engineering  viewpoint,  perhaps  the  over¬ 
arching  question  is  how  to  assess  the  dependency  of  aerodynamic  force  production  on  the  wealth  of 
parameters  -  Reynolds  number,  reduced  frequency,  wing  geometric  details  such  as  aspect  ratio  and  airfoil 
section,  the  role  of  structural  flexibility,  surface  roughness  and  so  forth.  One  is  overwhelmed  by  the  sheer 
range  of  independent  variables.  Then  there  is  the  question  of  how  much  fidelity  or  accuracy  is  required  in 
a  given  experiment  or  computation,  to  answer  the  basic  question  at  hand.  This  latter  question,  at  least, 
can  be  addressed  with  some  rigor  by  considering  what  we  call  “canonical  problems”. 


1.5  THE  NEED  FOR  CANONICAL  PROBLEMS 

By  “canonical  problems”  one  means  abstracted  test  cases  of  general  appeal,  with  clear  parameter  definitions, 
which  are  to  be  studied  by  several  computational  and  experimental  efforts,  by  different  research  groups 
working  towards  a  common  reporting  goal.  Objectives  include: 

1)  Cross-validation  of  the  various  methods; 

2)  Baselining  the  state  of  the  art; 

3)  Engendering  interest  in  the  subject  amongst  researchers  with  affinity  towards  this  area  but 
otherwise  lack  of  connection  to  the  specific  subject; 

4)  Establishment  of  a  known  test  case  for  archival  reference;  and 

5)  Providing  a  point  of  departure  for  individual  researchers  to  run  parameter  studies. 
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The  idea  is  that  once  a  common  baseline  is  achieved,  the  subsequent  parameter  studies  have  better 
contextual  motivation  and  are  on  firmer  ground  because  the  “correct”  answer  was  found  for  the  canonical 
case.  If,  on  the  other  hand,  significant  difference  is  found  in  results  for  the  canonical  case,  and  obvious 
errors  have  been  excluded,  then  the  canonical  problem  itself  requires  a  revisit  to  assess  flaws  in  state-of- 
the-art  computations  or  experiments,  before  further  work  can  be  justified. 

1.5.1  From  Fundamental  Research  to  Applications 

The  present  work  tends  toward  the  “fundamental”,  in  the  sense  that  no  flight  vehicle  or  prototype  in  nature 
is  analyzed,  and  the  methods  of  investigation  are  the  traditional  approach  in  laboratory  fluid  mechanics. 
However,  the  range  of  Reynolds  numbers,  motion  parameters  and  model  geometries  must  necessarily  be 
informed  by  the  applications  -  that  is,  MAV  configurations.  Pitching/plunging  airfoils  most  naturally  tend 
themselves  to  flapping-wing  MAV  applications,  for  example  in  problems  of  finding  good  thrust  efficiency 
due  to  leading  edge  suction. 

1.5.2  Role  of  Longitudinal  Problems 

Noting  the  complexity  of  flight  in  nature,  we  must  nevertheless  begin  with  suitably  simple  abstractions, 
ideally  those  stemming  from  classical  work  in  unsteady  aerodynamics  and  dynamic  stall,  but  suited  to  the 
Reynolds  numbers  and  reduced  frequencies  encountered  in  animal  flight  and  in  MAV  applications. 
A  logical  point  of  departure  is  longitudinal  problems,  where  the  airfoil  or  wing  is  oscillated  in  the 
longitudinal  plane;  that  is,  with  three  degrees  of  freedom:  an  in-plane  rotation,  or  pitch;  and  in-plane 
translation  normal  to  the  nominal  free-stream  direction,  or  plunge;  and  the  second  in-plane  translation, 
parallel  to  the  nominal  free-stream,  or  surge.  In  the  present  study,  interest  is  limited  to  pitch  and  plunge, 
but  as  discussed  below,  some  experimental  installations  have  a  slight  unintentional  but  unavoidable  surge. 


1.6  GENERAL  OUTLINE  OF  THE  REPORT 

We  consider  three  geometries:  a  low  Reynolds  number  airfoil  nominally  in  two  dimensions,  a  thin  flat  plate 
likewise  in  nominally  two  dimensions,  and  a  thin  flat  plate  of  aspect  ratio  2.  “Nominally  two  dimensions” 
means  for  experiments  that  the  model  spans  the  facility  test  section  from  wall-to-wall,  with  ideally  a  very 
small  gap  at  the  tips.  For  computations  the  meaning  is  more  ambiguous,  as  it  can  be  solving  the  Navier- 
Stokes  equations  (or  some  simplification)  in  a  plane,  or  in  some  volume  of  spanwise  extent  but  spanwise 
periodic  boundary  conditions. 

The  purpose  of  the  three  geometries  is  a  multi-way  comparison  of  sectional  geometry  effects  and  a 
preliminary  foray  into  3D  effects,  but  which  we  mean  whether  there  is  strong  departure  from  sectional 
aerodynamics  for  a  case  where  the  aspect  ratio  is  clearly  small. 

1.6.1  Listing  of  the  Cases  and  Philosophy  for  Why  They  Were  Chosen 

We  consider  the  following  kinematics  for  combined  pitch  and  plunge: 

Plunge:  h(t)  =  /z^  ccos(2;7- ft)  =  0.5ccos(p.5Uj  /  c) ;  and 

Pitch:  a{t)  =  -hA  cos(27r( ft  +  ^))  =  8°  +  8.42° cos(0.5f/^^ !  c  +  tt  12).  (3) 


1.6. 1.1  Airfoils 

Our  first  set  of  cases  features  the  Selig  SD7003  airfoil,  with  a  design  Reynolds  number  in  the  60000  - 
100000  range  [49].  This  was  also  motivated  by  the  steady-problem  studied  in  NATO  RTO  AVT-101, 
“MAV  Low  Reynolds  Number  Aerodynamics”  [50].  The  parameters  of  Eqn.  3  for  the  SD7003  airfoil  are 
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noted  schematically  in  Figure  1-8,  while  the  time  traces  of  effective  angle  of  attack  for  combined 
pitch-plunge  and  for  pure-plunge  are  given  in  Figure  1-9.  Our  choice  of  reduced  frequency, 
k  =  0.25  =coc/2Ua,  =  Trfc/Uao,  was  motivated  in  part  by  cruise-type  conditions  for  flapping  flight  of  birds. 
Although  the  Strouhal  number,  St  =  0.08,  is  below  the  range  for  maximum  propulsive  efficiency  for  most 
flyers  in  nature  [51],  the  present  flow  conditions  are  on  the  upper-end  of  the  dynamic-stall  literature  for 
helicopter  blade  applications  [5,8],  for  which  the  traditional  analytical  or  phenomenological  models  in 
aeronautics  tend  to  focus.  As  is  often  taken  in  applications  motivated  by  propulsive  efficiency  of  pitch/ 
plunge  [18],  pitch  leads  plunge  by  one  quarter  of  phase  (cp  =  0.25)  and  thus  the  airfoil  “feathers”,  with  the 
geometric  pitch  angle  partially  cancelling  the  plunge-induced  angle  of  attack,  arctan(^  /u^)- 

The  amplitude  of  pitch,  A  =  8.42°,  was  chosen  from  the  expression  ^  ^ _ X  is  the  ratio  of 

arctan(A„^Jt/„) 

pitch  angle  amplitude  to  the  peak  angle  of  attack  induced  by  the  plunge  motion;  we  chose  X  =  0.6,  which 
as  will  be  shown  below,  leads  to  shallow  dynamic  stall.  X  =  0,  on  the  other  hand,  is  a  pure-plunge,  which 
produces  a  strong  leading  edge  vortex,  and  is  more  akin  to  deep  dynamic  stall..  Variation  ofX  is  an  option 
for  parameter  studies  (not  pursued  here)  for  search  for  lift  and  thrust  efficiency,  while  keeping  Strouhal 
number  constant.  Alternatively,  Strouhal  number  can  be  varied  (by  changing  reduced  frequency  or 
reduced  amplitude)  and  X  varied  such  that  the  effective  angle  of  attack  history,  when  disregarding  pitch 
rate  effects,  is  kept  constant. 


Figure  1-8:  Schematic  of  SD7003  Airfoii  and  Parameters  of  Motion. 


h(t) 


Figure  1-9:  Motion  Kinematics  and  Effective  Angie  of  Attack  Time 
History  for  Pure-Plunge  and  Combined  Pitch-Piunge. 
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Figure  1-9  shows  the  plunge  trajectory  (green  curve),  total  effective  angle  of  attack  in  pure-plunge 
(black  curve),  pitch  geometric  angle  of  attack  (purple  curve)  and  total  effective  total  angle  of  attack  for 
combined  pitch-plunge  (blue  curves).  For  pitch-plunge,  the  total  effective  angle  of  attack  time-trace, 
ae,  straddles  the  static  stall  value  of  -  11°  [52];  this  is  just  the  sum  of  the  pitch  and  plunge  cosines  with 
appropriate  phase  shift.  But  ae,  can  be  taken  to  include  the  effect  of  pitch  rate,  which  depends  on  the  pitch 
pivot  point  location,  by  summing  all  of  the  components  inside  the  brackets  in  the  third  term  of  Equation  1 ; 
this  is  the  dashed  blue  curve  in  Figure  1-9.  The  difference,  vs.  disregarding  the  effect  of  pitch  pivot  point 
location  (solid  blue  curve)  is  a  phase  shift  of  -0.05  t/T.  With  inclusion  of  the  pivot  effect,  the  limits  on 
ae  become  2.03°  <  ae<  14.03°,  whereas  for  pure-plunge  they  are  -6.0°  <  ae<  22.0°. 

1.6.1.2  Flat  Plates  -  Wall-to-Wall 

The  second  set  of  cases  is  flat  plates,  nominally  of  3%  thickness  and  round  (circular)  edges.  In  some  sense 
the  flat  plate  is  a  more  fundamental  problem  than  the  airfoil,  because  apart  from  the  immediate  vicinity  of 
the  leading  edge,  there  is  no  pressure  gradient  due  to  the  model  geometry  itself  Plates  are  also  easy  to 
manufacture  for  experiments,  while  the  round  trailing  edge  somewhat  simplifles  the  computational  grid. 
And  if  the  plate  thickness  is  large  enough,  bending  and  other  structural  deflections  should  be  avoided. 

The  airfoil  cases  are  deeply  concerned  with  understanding  the  role  of  boundary  layer  transition  in  unsteady 
aerodynamics,  but  in  many  MAV  applications  transition  is  perhaps  of  secondary  importance.  Wings  tend  to 
be  thin,  with  sharp  edges;  and  transition  only  occurs  in  the  wake,  or  in  the  late  evolution  of  large  separated 
structures.  For  a  flat  plate  with  round  leading  edge  (as  opposed  to  say  a  “super-ellipse”),  presumably 
transition  would  occur  near  the  leading  edge,  and  would  be  fixed  across  a  broad  range  of  cases.  We  therefore 
repeated  the  pitch-plunge  and  pure-plunge  cases  for  the  airfoil,  but  now  with  the  plate.  But  as  with  the  airfoil, 
the  geometry  is  nominally  2D,  with  the  model  wall-to-wall  for  the  experimental  facilities,  and  2D  or 
spanwise  periodic  in  computations. 


1.6.1.3  Flat  Plates  of  Finite  Aspect  Ratio 

A  key  question  in  MAV  aerodynamics  is  the  role  of  aspect  ratio  effects  [3].  The  question  is  both 
elementary  and  subtle.  Of  course  one  is  interested  in  verifying  to  what  extent  lifting-line  or  lifting-surface 
theory  can  be  used  at  low  Reynolds  number  [53],  just  as  it  has  been  so  successful  for  manned  aircraft. 
Beyond  this,  there  is  the  question  of  how  the  tip  vortices  of  low  aspect  ratio  wings  might  interact  with  the 
LEV  to  maintain  flow  regularity  at  high  angle  of  attack  and  therefore  to  delay  stall  or  to  otherwise  provide 
lift  enhancement,  beyond  what  is  possible  in  2D  or  for  high  aspect  ratio  wings.  There  is  also  the  question 
of  how  lags  between  aerodynamic  forces  and  motion  kinematics  would  differ  in  3D  vs.  2D. 

Experimentally,  high  but  finite  aspect  ratio  wings  are  awkward  to  test.  Blockage  considerations  limit  span, 
while  high  aspect  ratio  means  very  small  chord,  which  in  turn  means  a  high  physical  oscillation  frequency 
for  a  given  reduced  frequency.  This  can  exceed  the  operating  parameters  of  motion  rigs,  especially  in  wind 
tunnels.  The  small  chord  is  also  problematic  for  optical  flowfield  velocimetry  methods.  In  computations, 
the  grid  would  have  to  be  large  in  the  spanwise  direction  to  capture  the  entire  wing  (or  half-wing)  and  the 
tip  region.  Thus  high-aspect  ratio  wings  are  more  difficult  to  study.  But  more  importantly,  we  deem  high 
but  finite  aspect  ratio  wings  to  be  not  very  interesting,  apart  from  very  applied  problems  such  as  the 
cruising  performance  of  large  birds.  Low  aspect  ratio  problems,  on  the  other  hand,  are  both  very 
interesting  and  geometrically  more  amenable  to  investigation  in  unsteady  conditions. 

We  therefore  somewhat  arbitrarily  chose  the  case  of  a  rectangular  wing  of  aspect  ratio  =  2,  with  the  same 
flat-plate  section  as  in  the  wall-to-wall  case,  and  round  edges  on  the  entire  perimeter.  Conveniently, 
for  AR  =  2  lifting-line  theory  and  slender  body  theory  give  the  same  inviscid  lift  curve  slope,  tt. 
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Chapter  2  -  FACILITIES  AND  METHODS 

2.1  GENERAL  COMMENTS 

Here  we  describe  the  various  experimental  facilities  and  computational  methods  used  by  the  Task  Group. 
As  the  entire  approach  is  predicated  on  cross-validation  between  the  various  facilities,  between  the  various 
codes,  and  between  the  two,  it  is  important  to  recite  the  capabilities  and  limitations  of  each.  For  example, 
a  given  facility  might  have  higher  or  lower  turbulence  intensity  than  another,  or  a  given  model  motion  rig 
might  have  higher  or  lower  fidelity  of  motion  (vibration,  elastic  deflection,  motion  overshoot,  etc.) 
than  another.  These  may  result  in  transition  effectively  further  upstream.  Further,  blockage  may  differ, 
as  would  the  role  of  “3D  effects”,  support  interference,  flow  non-uniformity,  and  on  and  on.  Similarly, 
the  computational  approaches  differ  in  resolution  and  mesh  type,  2D  vs.  3D,  treatment  of  boundary 
conditions,  modelling  of  the  far-field,  type  of  filtering  and  averaging,  turbulence  and  transition  modelling, 
and  so  forth.  While  it  is  impossible  to  rigorously  note  all  of  these  features,  our  aim  is  to  be  sufficiently 
descriptive  to  enable  placement  in  context  of  the  various  results,  such  that,  for  example,  a  large  laminar 
separation  in  a  “quiet”  tunnel  and  a  much  smaller  separation  in  a  “noisier”  tunnel  for  nominally  identical 
conditions  could  be  ascribed  to  the  right  interpretation  of  the  respective  tunnel  properties.  Towards  that  end, 
the  more  experiments  and  computations  on  identically-posed  runs,  the  better. 

We  note  separately  the  various  facilities  using  water  as  the  working  fluid,  the  facilities  using  air,  the  force 
measurement  schemes,  and  the  computational  methods. 


2.2  WATER  TUNNELS  AND  WATER  TOW  TANKS 

For  unsteady  aerodynamics,  the  advantage  of  liquid-based  facilities,  as  opposed  to  wind  tunnels,  is  that  for 
the  same  reduced  frequency  the  physical  frequency  of  motion  is  much  lower.  This  reduces  loads  on 
motion  mechanisms  and  reduces  the  data  acquisition  rate  requirements  for  diagnostics.  Disadvantages 
include  corrosion  and  conductivity  of  water,  which  complicate  choices  for  materials  and  for  force 
balances.  Water  tunnels  also  tend  to  have  higher  turbulence  intensity  than  do  the  comparable  wind  tunnels. 
In  theory  tow  tanks  are  the  best  choice,  with  zero  turbulence  intensity,  but  have  the  disadvantage  of 
requiring  drive  mechanisms  and  long  settling  times  between  successive  tests. 

2.2.1  University  of  Michigan  Water  Tunnel 

The  University  of  Michigan  (UM)  water  tunnel  has  a  test  cross-section  61  cm  wide  by  61  cm  high  and  free 
stream  speed  range  of  5  cm/s  -  40  cm/s,  with  a  turbulence  intensity  of  approximately  1%,  as  measured  by 
taking  free-stream  PIV  images  over  a  time -window  of  hours.  However,  the  short  time  turbulence  intensity 
computed  by  averaging  over  a  time  period  on  the  order  of  the  test  section  residence  time  is  less  than  0.5%. 
The  larger  value  of  the  long  time  average  is  associated  with  very  low  frequency  sloshing. 

SD7003  and  flat  plate  airfoil  models  with  152  mm  chord  were  tested.  The  SD7003  airfoil  model  was 
fabricated  with  stereo  lithography  and  a  transparent  resin  (DSM  Somos  11122)  to  minimize  laser  reflection 
at  the  model  surface,  while  the  flat  plate  model  was  built  of  stainless  steel,  with  a  polished  surface  to  reduce 
laser  glare.  Both  models  were  hung  vertically,  spanning  the  depth  of  the  test  section  and  terminating  at  1  mm 
from  the  bottom,  with  an  end  plate  at  the  free-surface.  This  cantilevered  mounting  scheme  results  in  a 
maximum  tip  deflection  of  -0.5  cm  for  the  airfoil,  and  0.1  cm  for  the  plate.  Pitch  is  produced  by  a  rotary 
stage  (Velmex  B4872TS  Rotary  Table).  One  linear  traverse  (Velmex  20-inch  BiSlide)  produces  plunge, 
and  a  second  (Velmex  40-inch  BiSlide)  produces  streamwise  motion  or  surge  (not  used  here).  All  are 
mounted  above  the  test  section  and  controlled  by  a  Velmex  VXM-1-1  motor  controller.  The  experimental 
setup  is  illustrated  in  Figure  2-1. 


RTO-TR-AVT-149 


2-1 


FACILITIES  AND  METHODS 


ORGAISriZATION 


Nd-YAG 


Optics 


Bi-slides 

Rotary  stage 


Pitch  motion 


Top  View 


A 

Plunge 

motion 

V 


Water  surface  -  -  -  . 


Side  View 


Freesiream 

Velocity 


<■ 


Rotary  stage 


Figure  2-1:  University  of  Michigan  Water  Tunnei:  Motion  Rig  Schematic  (ieft) 
and  instaiiation  Above  Tunnei  Test  Section  (right). 


Measurements  consisted  of  dye  flow  visualization  and  particle  image  velocimetry  (PIV).  The  dye  injection 
system  is  a  probe  rake  with  seven  evenly  spaced  dye  streams  introduced  half  chord  in  front  of  the  leading 
edge  of  the  model.  Two  syringe  pumps  are  used  to  match  the  dye  speed  at  the  injection  point  and  the  water 
channel  flow  speed  to  minimize  disturbances.  The  dye  streams  are  held  fixed  and  the  airfoil  through  the 
dye  streams.  The  resulting  streaklines  capture  flow  direction  and  recirculation  zones. 

The  PIV  system  includes  a  double-pulsed  Nd:YAG  laser  (Spectra  Physics  PIV  300)  and  two  dual  frame 
digital  cameras  (Cooke  Corp.  PC0.4000).  PIV  seeding  was  with  3  pm  diameter  Titanium  Dioxide  particles. 
A  small  amount  (8  drops)  of  a  dispersant  (DARVAN  C-N,  Vanderbilt)  was  added  to  the  5,000  gallon  water 
tunnel  to  improve  particle  distribution  and  suspension,  and  to  avoid  clumping.  The  cameras  were  installed 
underneath  the  test  section  and  were  equipped  with  Nikon  105-mm  Micro-Nikkor  lenses  to  produce  a 
magnification  of  approximately  25  pixels/mm  and  a  field  of  view  of  160  by  107  mm.  Time  between 
exposures  was  adjusted  to  produce  a  nominal  particle  displacement  of  8  pixels  at  the  free  stream  velocity  for 
all  cases.  To  capture  the  h  =  0.5  c  plunge  motion  and  to  avoid  shadowing  of  the  field  of  view,  the  leading 
edge  and  trailing  edge  of  the  model  were  imaged  separately,  with  an  overlap  region,  and  stitched  together 
before  processing.  The  accuracy  of  the  axial  Velmex  BiSlide  traverse  is  0.00635  mm,  giving  approximately 
1/6  of  a  pixel  of  image  overlap  uncertainty. 

200  PIV  pairs  were  taken  at  each  phase,  with  the  Nd:YAG  laser,  CCD  cameras,  rotary  stage,  and  BiSlide 
were  phase-locked.  Every  run  recorded  12  periods  and  only  the  last  10  used  to  compute  the  phase  averages. 
Recording  was  initiated  by  the  PIV  system  data  acquisition,  which  triggered  the  motion  controller.  The  PIV 
laser  timing  and  the  motion  period  were  matched  with  an  accuracy  of  0.1  ms  for  a  typical  period  of 
approximately  10  s.  This  produced  a  slight  discrepancy  in  the  model  position  between  the  first  image  and  the 
last  image  at  phases  with  large  plunge  speed  (7  pixels  max,  or  0.28  mm).  This  discrepancy  resulted  in  an 
elimination  of  a  datum  point  near  the  model  surface. 

PIV  images  were  analyzed  with  in-house  MATLAB-based  PIV  software.  Particle  displacement  is 
determined  in  two  passes  using  cross-correlation  of  displaced  interrogation  windows  and  Gaussian  fit  for 
sub-pixel  fit.  The  first  pass  uses  64  x  64  pixel  windows  with  uniform  8  pixel  window  displacement. 
The  second  pass  uses  32  x  32  pixel  windows  with  window  displacement  from  the  PIV  calculation  in  the 
first  pass,  giving  an  approximate  spatial  resolution  of  0.64  mm. 
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Several  validation  criteria  were  applied  to  the  measured  particle  displacements.  The  displacement  must  be 
within  a  predetermined  range  of  values  in  the  x-  and  y-directions.  The  range  of  values  in  the  first  pass  is 
fairly  large,  to  capture  the  large  range  of  particle  displacements  found  near  the  model  surface;  and  small 
(10  pixels  displacement)  in  the  second  pass.  A  median  filter  is  used  to  find  the  particle  displacement  at  the 
points  where  the  PIV  validation  failed,  and  to  remove  outliers.  A  square  grid  with  16  pixel  spacing  was 
used  for  all  the  images.  Near  the  surface  of  the  airfoil,  data  points  within  32  pixels  from  the  boundary  were 
discarded  because  the  interrogation  window  would  include  pixels  inside  the  model.  This  corresponds  to 
two  data  points  in  the  measurement  grid.  To  further  remove  outliers  based  on  large  sample-to-sample 
fluctuations,  a  3-sigma  filter  was  used.  First,  the  sample  mean  and  standard  deviation  from  the  200  PIV 
images  at  each  phase  were  calculated  at  each  point  in  the  flowfield.  Then  with  the  sample  mean  and 
standard  deviation  value,  each  datum  point  was  revisited  and  the  value  was  compared  to  ±3  standard 
deviation  of  the  mean  value.  The  datum  point  was  discarded  if  it  lied  outside  the  3 -sigma  region  and  the 
remaining  data  were  used  to  calculate  a  more  accurate  mean  flow  and  turbulence  statistics.  Out  of 
200  images  recorded,  approximately  195  images  were  used  for  phase  averages  for  relatively  attached 
flowfields.  For  flows  with  large  separation  regions,  approximately  180  images  were  retained  after  the 
3 -sigma  filter.  The  calculated  95%  confidence  interval  for  the  mean  and  turbulence  statistics  is  tabulated 
in  Table  2-1. 


Table  2-1 :  Summary  of  95%  Confidence  Intervals  for  UM  PIV  Data 
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±  0.07% 

-v  ±  0.7% 

±  0.13% 

±  1.25% 

In  Table  2-1  the  values  of  the  95%  confidence  intervals  are  normalized.  Velocities  such  as  u  are  normalized 
by  the  free  stream  value,  and  statistics  such  as  uv'  are  normalized  by  the  free  stream  velocity  squared. 
For  each  model  and  Reynolds  number  combination  two  values  of  the  95%  confidence  interval  are  reported; 
the  mean  and  the  maximum.  The  mean  is  the  mean  of  the  95%  confidence  intervals  computed  for  all  the 
phases  in  a  particular  case.  The  maximum  simply  reports  the  maximum  95%  confidence  interval  in  the  entire 
flowfield  for  all  phases.  This  value  indicates  the  worst  possible  uncertainty  on  the  calculated  mean  value  for 
a  given  airfoil  kinematics.  In  general,  the  mean  and  the  maximum  confidence  intervals  at  motion  phases 
where  the  flow  is  largely  attached  were  an  order  of  magnitude  less  than  at  phases  with  large  separation. 

2.2.2  U.S.  Air  Force  Research  Laboratory  Water  Tunnel 

The  U.S.  Air  Force  Research  Laboratory’s  horizontal  free-surface  water  tunnel  has  4:1  contraction  and 
46  cm  wide  by  61  cm  high  test  section.  Free-stream  speed  range  is  3  -  45  cm/s.  A  photograph  of  the  tunnel, 
schematic  of  the  pitch-plunge  rig  and  installation  of  the  SD7003  airfoil  model  in  the  tunnel  are  shown  in 
Figure  2-2.  In  the  photograph,  the  model  is  inside  the  test  section,  but  the  test  section  glass  sidewalls  are  not 
visible. 
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Figure  2-2:  AFRL  Horizontal  Free-Surface  Water  Tunnel  (left);  Schematic  of  Pitch-Plunge  Rig  and 
Airfoil  Model  (middle);  SD7003  Airfoil  Installed  Inside  Test  Section,  Showing  Smooth  Suction- 
Side  of  Airfoil  (top  right)  and  Plunge  Rod  Coupling  on  Pressure-Side  of  Airfoil  (bottom  right). 
Black  arrow  in  bottom-right  points  to  dye  injection  exit  port  location. 


Water  tunnel  turbulence  intensity  (u  and  v  components)  as  determined  from  PIV  free-stream  data  and 
confirmed  by  laser-Doppler  velocimetry,  is  estimated  at  0.4%  at  C/oo  =  30  -  40  cm/s.  A  surface  skimmer 
plate  at  the  entrance  to  the  test  section  and  a  sealed  lid  of  the  intake  plenum  (visible  as  plywood  cover  in 
Figure  2-2)  damp  sloshing  in  the  tunnel,  that  could  otherwise  have  been  be  excited  as  a  first-mode  shallow 
water  wave,  and  would  have  resulted  in  low-frequency  (~0.16  Hz)  sinusoidal  variation  in  f/oo.  More  detail 
on  the  facility  is  given  in  01  et  al.  [52]. 

The  tunnel  is  fitted  with  a  3 -component  oscillator  rig,  with  three  electric  linear  motors,  controlled  through 
a  Gain  DMC  4040  Ethernet  controller.  Pitch  and  plunge  are  made  possible  by  a  pair  of  motors  mounted 
vertically  on  a  plate  above  the  tunnel  test  section.  Each  motor  actuates  a  vertical  “plunge  rod”,  which 
connects  via  a  bushing  to  the  model  at  a  fixed  pivot  point  on  the  model  chord,  in  the  test  section  vertical 
plane  of  symmetry.  The  upstream  plunge  rod  is  constrained  to  move  purely  vertically,  whereas  the 
downstream  plunge  rod  is  allowed  to  pivot  in  the  test  section  vertical  plane  of  symmetry.  The  desired 
angle  of  attack  and  vertical  position  time  history  of  the  model  are  converted  to  position  commands  for 
each  linear  motor.  This  allows  for  single  degree-of- freedom  motions  such  as  pure -pitch  about  a  prescribed 
fixed  pivot  point,  or  pure-plunge.  Pitch  and  plunge  can  be  combined,  and  the  pitch  pivot  point  can  be 
varied  by  suitable  choice  of  phase  and  amplitude  difference  in  trajectory  of  front  or  rear  plunge  rod.  For  all 
cases  where  the  pitch  pivot  point  is  not  coincident  with  the  bushed  end  of  the  front  plunge  rod,  there  will 
be  a  small  parasitic  streamwise  displacement  of  the  model,  which  would  be  unavoidable  unless  the  front 
plunge  rod  were  to  be  allowed  to  pivot  similarly  to  the  downstream  one.  A  third  degree  of  freedom,  surge, 
is  made  possible  by  a  larger  linear  motor  mounted  horizontally  aft  of  the  pitch-plunge  carriage.  It  is  used 
to  remove  the  aforementioned  parasitic  streamwise  displacement,  as  was  necessary  for  aft-mount  sting 
experiments  with  an  AR  =  2  plate  described  in  Chapter  6. 

Particle  image  velocimetry  measurements  were  taken  with  a  PCO  4000  1  IMpix  camera  run  in  single  A-D 
mode,  triggered  from  an  external  pulse  train  derived  from  the  position  encoder  of  the  motion  rig, 
thus  allowing  for  selection  of  motion  phase  at  which  to  acquire  data,  and  for  phase  averaging.  Velocity 
data  were  acquired  once  per  period  of  oscillation,  for  a  sequence  of  120  periods.  The  first  5  periods  were 
removed  from  each  data  set,  and  the  remaining  115  velocity  data  sets  were  ensemble  averaged  for  each 
phase  of  motion.  Vorticity  was  calculated  by  explicit  differentiation  of  cubic  spline  fits  to  the  velocity 
field  (as  discussed  by  Willert  and  Gharib  [54]).  PIV  resolution  was  typically  156  pixels/cm,  which  for 
32  X  32-pixel  windows  with  16  x  16  overlap  gives  148  vectors  per  chord  length.  Because  of  laser 
reflections  from  the  model  surface  (polished  stainless  steel,  spray-painted  flat-black)  and  lack  of 
corrections  for  PIV  windows  which  at  least  partially  intersect  with  the  model  surface,  data  in  the  first  row 
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of  PIV  windows  adjacent  to  the  airfoil  surface  are  not  reliable.  Light  sheet  thickness  was  approximately 
1.5  mm  max,  though  the  large  field  of  coverage  (up  to  45  cm  across  in  some  cases)  makes  precise 
collimation  of  the  light  sheet  difficult.  To  obtain  the  “most  2-dimensional”  flowfield,  the  PIV  light  sheet 
was  placed  at  the  3/4  span  location;  that  is,  approximately  halfway  between  the  plunge  rods  and  the  tunnel 
wall,  and  0.05  c  inboard  of  the  dye  injection  slot.  More  detail  on  the  rig  operation  and  discussion  of 
experimental  error  is  given  in  01  et  al.  [55]. 

For  dye  injection,  a  wand  with  0.5  mm  internal  diameter,  injecting  concentrated  blue  food  colouring, 
was  fitted  inside  the  model  and  exits  flush  with  the  outer  mould  lines  near  the  airfoil  leading  edge,  firing 
approximately  wall-normal.  The  dye  exit  location  is  not  visible  in  Figure  2-2,  but  is  marked  by  the  black 
arrowhead  in  the  bottom-right-hand  portion  of  Figure  2-2. 

2.2.3  Cambridge  University  Towing  Tank 

The  Cambridge  University  Engineering  Department  (CUED)  towing  tank  is  7  m  long  and  has  a  1  m  x  1  m 
cross-section  with  a  central  2  m  long  Perspex  test  section.  A  three-dimensional  model  and  photo  of  the 
waving  arm  mechanism  is  shown  in  Figure  2-3,  where  the  entire  mechanism  is  shown  suspended  from  the 
tow  tank  carriage  on  four  vertical  struts.  The  Perspex  sides  and  the  water  volume  of  the  tank  are  shown  to 
provide  a  sense  of  scale.  The  water  depth  is  0.8  m. 


Figure  2-3:  Schematic  of  Section  of  CUED  Water  Tow  Tank  and  “Waving  Arm”  Setup. 

The  waving  wing  mechanism  was  designed  to  mimic  the  translational  phases  of  an  insect  wing  stroke, 
extending  the  propeller  experiments  to  include  wing  starting  and  stopping.  The  wing  is  started  from  rest 
and  swept  through  a  wing  stroke  angle  0  up  to  90  degrees,  with  or  without  a  “freestream”  generated  by 
translating  the  mechanism  via  a  towing  tank  carriage.  The  wing’s  angular  velocity  can  be  programmed  to 
any  profile. 
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The  mechanism  is  mounted  on  a  stainless  steel  plate  suspended  just  below  the  water  surface.  The  waving 
motion  is  controlled  by  a  servomotor  and  gearbox  programmable  through  Lab  View,  and  a  slotted  optical 
switch  is  mounted  near  one  end  of  the  axle  to  confirm  wing  position.  A  force  balance  can  be  mounted  on 
the  wing  to  provide  unsteady  lift  and  drag  force  measurements.  The  distance  from  the  wing  root  (defined 
as  the  bottom  of  the  skim  plate  when  the  wing  is  vertical)  to  the  axis  of  rotation  is  10%  of  the  span. 
The  angle  of  attack  is  selectable  from  0  to  45  degrees  in  5  deg  increments.  The  wing  is  a  2.5%  thick 
carbon  fibre  flat  plate  with  rounded  edges,  a  chord  of  0.125  m,  and  an  aspect  ratio  of  4  as  measured  from 
the  free  surface  to  wingtip.  A  90  deg  wing  stroke  is  equivalent  to  5.37  chords  of  travel  at  3/4  span.  At  a 
15  deg  angle  of  attack  blockage  of  the  towing  tank  is  2%. 

The  wing  stroke  was  programmed  using  three  different  velocity  profiles:  linear,  sinusoidal,  and  exponential. 
The  wing  was  accelerated  such  that  it  reached  its  maximum  velocity  after  0.10  (exponential  only),  0.25, 
or  0.60  chord-lengths  of  travel  for  each  profile.  Velocity  profiles  were  symmetric  such  that  the  wing 
decelerated  in  the  same  way  before  reaching  a  maximum  stroke  angle  of  90  degrees.  The  maximum 
velocity  was  chosen  as  that  which  gives  a  local  wing  Reynolds  number  of  60000  at  3/4  span. 

2.2.4  ONERA  Tow  Tank 

The  ONERA  (Lille)  water  towing  tank  has  a  carriage  system  with  towing  speed  of  between  0.1  m/s  and 

1.4  m/s.  Speeds  below  0.1  m/s  are  unstable,  while  at  speed  above  1.4  m/s  the  available  test  duration  is  too 
short.  The  towing  tank  can  alternatively  be  run  as  a  water  tunnel,  but  the  flow  is  not  uniform  and  is 
unsteady  with  a  large  turbulence  level  (>3%).  So  the  results  presented  here  are  only  obtained  in  towing- 
tank  mode.  Tests  with  hot  wires  in  the  water  show  that  the  residual  velocities  before  test  start  are  too  low 
to  be  measured  (<1  mm/s).  As  the  facility  was  built  quite  some  time  ago,  the  servo  command  of  the 
carriage  velocity  is  not  as  efficient  as  desired.  The  length  needed  to  reach  operating  velocity  is  a  couple  of 
periods.  As  the  unsteady  hydrodynamic  effects  are  important  during  these  first  periods,  the  measured 
forces  on  the  model  are  retained  only  are  after  a  couple  of  periods.  To  ensure  that  the  frequency  of 
movement  is  constant  and  to  phase  the  measurement  at  the  same  position  from  each  period,  an  optical 
trigger  system  is  used  on  one  of  the  rotating  wheels. 

The  model  oscillation  mechanism  consists  of  two  crank  wheel  systems  with  a  phase  lag  between  the  two 
wheels.  The  amplitude  in  pitch  depends  on  the  difference  of  phase  between  the  wheels.  The  mechanisms 
were  conceived  for  flap  and  pitch  motions  as  well  as  plunge  and  pitch  motions,  and  therefore  require 
adjustment  to  produce  true  pitch-plunge.  Approximation  to  true  trigonometric  motion  depends  on  the  crank 
stroke  to  connecting  rod  length  ratio,  which  is  large  -  and  hence  the  approximation  is  close.  To  reduce 
vibrations,  a  belt  drive  system  is  used.  Servo  control  motors  generate  excitation  of  the  mechanism  at 
frequencies  higher  than  the  periodic  oscillations;  specific  transfer  functions  of  the  servo  mechanism  are 
necessary  to  eliminate  excitations  at  a  frequency  near  the  structural  resonance.  As  the  system,  with  a  constant 
rotation  rate  motor,  is  not  controlled  in  a  closed  loop,  the  movement  is  forced  by  the  mechanism. 
The  mechanism  and  model  installation  are  shown  in  Figure  2-4. 


2-6 


RTO-TR-AVT-149 


NATO 

OTAN 


FACILITIES  AND  METHODS 


Figure  2-4:  ONERA  Tow  Tank  SD7003  Wing  Instailation;  View  of 
Model  (left)  and  Oscillation  Rig  Atop  the  Towing  Tank  (right). 


To  avoid  free  surface  effects,  that  are  important  at  relatively  high  velocity  (>0.5  m/s),  a  flat  plate  moving 
with  the  wing  can  be  fixed  on  the  strut  supporting  the  balance.  For  quasi-2D  runs  the  tip  of  the  wing  can 
be  close  to  the  tank  ceiling,  but  because  the  tank  floor  is  somewhat  uneven,  the  tip  gap  can  not  be  too 
small.  The  physical  aspect  ratio  of  the  SD7003  airfoils  tested  in  the  tank  was  7  or  14. 


2.3  WIND  TUNNELS 

2.3.1  Technical  University  of  Braunschweig  Wind  Tunnel 

The  TUBS  facility  is  an  atmospheric  Eiffel-type  wind  tunnel  with  a  closed  test  section  (Figure  2-5). 
The  settling  chamber  with  nozzle,  the  diffuser  and  the  motor  mounting  are  made  from  glass  fibre 
reinforced  epoxy  and  the  test  section  is  transparent.  A  3  kW  motor  drives  the  tunnel  in  suction  mode. 
The  wind  tunnel  is  installed  in  a  room  which  allows  circulation  of  the  flow.  Motor  and  fan  are  decoupled 
from  the  rest  of  the  tunnel. 
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Figure  2-5:  Low-Noise  Wind  Tunnei  of  Technische  Universitat  Braunschweig. 


The  wind  tunnel  nozzle  has  a  rectangular  cross-section  with  a  Boerger-type  axial  area  distribution. 
The  contraction  ratio  is  about  16  to  achieve  low-turbulence  flow  in  the  test  section.  The  exterior  wall  is 
covered  with  acoustic  foam  plates  to  suppress  noises  reflections.  At  the  inlet  to  the  settling  chamber 
several  devices  are  installed  to  reduce  turbulence  and  to  achieve  a  uniform  flow.  At  first  the  air  flows 
through  a  30  mm  thick  fleece  mat.  It  then  passes  through  a  honeycomb  flow  straightener  and  finally 
through  a  fine  woven  screen  to  decrease  the  size  of  the  turbulent  structures.  The  test  section  has  a  height  of 
600  mm,  a  width  of  400  mm  and  a  length  of  1500  mm.  A  vane  anemometer  is  used  to  control  the  speed  of 
the  flow  in  the  test  section.  The  measured  velocity  variations  in  the  core  flow  are  below  1%. 
The  turbulence  level  was  measured  by  means  of  the  hot  wire  anemometry.  The  turbulence  level  measured 
along  the  vertical  symmetry  axis  of  the  test  section  is  around  0.16%  at  10  m/s. 

The  2D  airfoil  motion  apparatus  of  TU  BS  is  mounted  around  the  test  section.  The  task  of  the  apparatus  is 
to  simultaneously  perform  plunging  and  pitching  motion:  A  3  kW  motor  drives  a  4:1  ratio  transmission 
which  transfers  the  rotation  motion  to  two  interdependent  systems:  The  plunging  motion  system  and  the 
pitching  motion  system.  Plunging  motion  is  created  by  a  sliding  carriage  connected  to  a  balanced 
flywheel.  The  sliding  carriage  has  a  lightweight  design.  A  timing  belt  drives  a  second  sliding  carriage  for 
supplying  the  pitching  motion.  The  apparatus  allows  adjusting  plunge  amplitudes  of  100  mm,  pitch  angles 
up  to  25  deg,  and  the  motion  frequency  may  be  as  large  as  5  Hz. 

Boundary  layer  resolution  for  phase-locked  PIV  measurements  is  obtained  from  1 1  measurement  windows 
along  the  upper  surface  of  SD  7003  airfoil,  according  to  Figure  2-6.  For  each  of  the  selected  phase  angles 
at  least  800  image  pairs  are  obtained.  Lack  of  a  direct  camera  view  of  the  PIV  light  sheet  in  the  streamwise 
plane  resulted  in  the  need  for  askance  viewing,  which  fits  naturally  with  rotational-type  stereoscopic  PIV 
with  Scheimpflug  camera-lens  adapters  to  remove  astigmatism. 


2-8 


RTO-TR-AVT-149 


FACILITIES  AND  METHODS 


0.1 

0.05 

0 

a 

-0.05 

-0.1 


v«iciaty_Llinf 


1.55 


HO 


125 


1,10 


0.95 


0.60 


0.65 


0.50 


0.55 


ozo 


0.05 


Flapping  Flight 
Motion  Apparatus 


.  .  .  . .  ^  Double  pulsed 


Nd:YAG  laser 


Motor  control 


PIV  aquisition 


Figure  2-6:  Arrangement  of  PIV  Interrogation  Windows  (left)  and  Stereo 
PIV  Setup  Used  for  the  Flowfield  Measurements  of  TUBS  (right). 


Accurate  calibration  is  critical  to  successful  stereo  PIV.  A  calibration  plate  is  therefore  positioned  coplanar 
to  the  laser  light  sheet,  and  set  a  several  displacements  normal  to  the  light  sheet.  The  images  will  not  be 
sharp  for  every  position  y  since  the  cameras  are  only  focused  on  one  certain  plane  and  so  the  aperture  of 
the  camera  objective  has  to  be  optimized.  The  calibration  concerns  capturing  images  of  the  grid  on  the 
calibration  plate.  Note  that  one  optical  path  from  one  camera  has  to  pass  the  plate.  When  performing  PIV 
measurements,  i.e.  capturing  particle  images  of  the  laser  light  sheet,  there  is  no  plate  in  the  optical  path, 
and  the  optical  path  is  slightly  different.  This  difference  is  called  disparity.  It  is  estimated  by  correlating  a 
single  particle  image  seen  by  the  different  views  of  the  two  cameras.  The  disparity  map  is  used  afterwards 
to  correct  the  calibration  matrix  by  analytical  formulations. 

PIV  data  processing  for  a  moving  model  requires  several  steps  in  addition  to  the  usual  PIV  cross¬ 
correlations.  One  issue  is  “wobbling”  of  the  airfoil  surface,  because  the  oscillatory  motion  is  not  exact  from 
period  to  period.  A  wobble  correction  is  therefore  performed.  When  the  laser  light  sheet  reaches  the  airfoil 
surface,  the  reflection  line  is  visible  on  the  camera  images.  Due  to  the  phase  locked  imaging,  this  reflection 
line  which  indicates  the  position  of  the  airfoil  should  be  always  at  the  same  location.  However,  the  reflection 
line  is  wobbling  about  0.5  millimetres  in  the  camera  images.  This  wobbling  has  to  be  removed  since  the 
airfoil  has  to  be  at  the  same  position  for  the  later  ensemble  averaging  procedure  of  the  vector  fields. 
The  required  shifting  is  performed  with  a  correlation  algorithm. 

In  the  next  step  of  improving  the  particle  image  quality,  all  particle  image  pairs  are  averaged.  Hence, 
the  particle  information  vanishes,  but  the  laser  reflection  line  and  the  average  noise  in  the  image 
background  remain.  This  average  image  is  now  subtracted  from  each  particle  image.  Consequently, 
the  intensity  of  the  reflection  line  and  the  background  noise  decreases,  there  is  a  better  signal  to  noise 
ratio.  Furthermore,  the  areas  in  the  particle  images  inside  the  airfoil  without  information  are  masked  out. 

Having  completed  several  image  pre-processing  techniques  to  improve  the  particle  image  quality, 
the  particle  displacement  evaluation  is  executed  in  the  next  step  using  a  cross-correlation  scheme.  A  multi¬ 
pass  interrogation  scheme  is  applied  with  decreasing  interrogation  window  size  (from  128  x  128  pixels  down 
to  32  X  32  pixels),  50%  overlap  and  elliptical  weighting  function.  The  resulting  set  of  vector  fields  for  each 
measurement  window  is  later  post  processed.  This  is  necessary  to  filter  out  non-physical  vectors  in  the  vector 
fields  which  would  corrupt  the  results  of  the  ensemble  averaging  procedure.  For  this  purpose  a  correlation 
peak  ratio  filter  is  employed  with  a  peak  ratio  of  3.  Furthermore,  iterative  remove  and  replace  filtering  is 
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used.  The  vector  of  a  certain  interrogation  window  was  only  considered  for  ensemble  averaging  if  at  least 
800  image  pairs  produced  valid  vector  data.  More  details  are  found  in  Bansmer  et  al.  [56,57]. 

2.3.2  DLR  Wind  Tunnel  PIV  Measurements 

The  German  Aerospace  Centre  (DLR),  Institute  for  Aerodynamics  and  Flow  Technology,  pursued  a 
combination  of  particle  image  velocimetry  and  surface  photogrammetry  for  a  range  of  flexible  wing  pitch- 
plunge  experiments  in  the  TUBS  wind  tunnel,  already  described  above.  The  surface  photogrammetry 
enabled  determination  of  deflection  time  history  of  flexible  structures  and  thus  fluid-structure  interaction, 
but  for  the  nominally  rigid  wings  of  the  present  study,  this  capability  was  not  of  direct  priority. 
Accomplishing  both  measurements  simultaneously  requires  care  in  surface  coatings  and  optical  filters  of 
the  imaging  cameras;  unfortunately  this  sometimes  leads  to  camera  pixel  saturation  from  laser  reflections, 
exacerbating  the  usual  problem  of  an  “optical  boundary  layer”  plaguing  PIV  at  the  model  surface. 

The  same  optical  setup  is  used  for  measuring  the  model  position  and  orientation  and  to  measure  the  wing 
shape  and  deformation  (Figure  2-7).  Two  JAI  Al  cameras  (1392  x  1024  pix)  are  arranged  above  the  top 
window  of  the  test  section  symmetrically  along  the  centre  line  of  the  model.  The  lenses  (f  =  16  mm)  are 
chosen  such  that  the  model  is  always  completely  captured  from  both  cameras  at  viewing  angles  20°. 
The  upper  model  surface  is  homogenously  illuminated  from  above  by  an  extended  beam  from  a  25  mJ 
Nd:YAG  pulse  laser.  The  PIV  setup  has  two  laser  light  sheets  to  capture  the  flowfleld  in  a  streamwise 
plane  below  and  above  the  moving  wing  at  a  fixed  span  position  of  b/4.  Two  double  pulsed  Nd:YAG  laser 
heads  (Big  Skye  CFR-400)  each  providing  about  150  mJ  light  energy  per  pulse  are  located  behind  the  top 
and  bottom  window  of  the  test  section.  Two  separate  PCO-1600  double  shutter  cameras  (1600  x  1200  pix) 
are  used  to  observe  the  flowfleld  below  and  above  the  wing  at  a  right  angle  to  the  light  sheets.  The  optical 
measurements  were  performed  simultaneously  at  10  different  phase  positions  of  the  sinusoidal  movement 
of  the  model.  For  each  run  1050  images  were  taken  by  each  camera  allowing  a  phase  averaging  of  the 
results  using  105  samples.  For  data  processing,  the  “Point  Tracker”  and  “Strain  Master”  of  the  DaVis  PIV 
software  are  used;  more  details  are  given  in  Kirmse  and  Wagner  [58],  and  Konrath  et  al.  [59]. 
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IPCT-cameras  iPCT-illuinfnation  laser 


2.3.3  TU  Darmstadt  Wind  Tunnel 

The  TU  Darmstadt  experimental  setup  consists  of  a  pair  of  pitch-plunge  motion  rigs  and  an  Eiffel-type  wind 
tunnel  with  test  section  of  cross-section  of  45  cm  by  45  cm  and  a  length  of  2  m.  The  tunnel  has  a  contraction 
ratio  of  24:1  with  five  turbulence  filters  in  the  settling  chamber  and  produces  turbulence  levels  on  the  order 
of  1.0%  at  the  test  speeds  of  interest  here.  A  schematic  of  the  wind  tunnel  is  shown  in  Figure  2-8.  The  free- 
stream  velocity  was  controlled  via  closed-loop  control,  with  the  tunnel  speed  input  obtained  from  a  hot-wire 
anemometer  (Dantec  Dynamics  A/S  type  55P11)  positioned  at  the  entrance  of  the  test  section.  The  hot-wire 
anemometer  was  calibrated  for  each  new  set  of  measurements  using  a  miniature  vane  anemometer. 
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Figure  2-8:  Schematic  of  Wind  Tunnei;  Note  intake  (A-A)  and  Test  Section  (C-C)  Cross-Sections. 


A  commercial  PIV  system  (Dantec  Dynamics  A/S)  consisted  of  a  Nd:YAG  {X  =  532  nm)  Litron  dual-cavity 
laser  with  a  maximum  power  output  of  135  mJ  per  cavity  and  two  10-bit  FlowSense  2M  CCD  cameras  each 
with  a  1600  x  1200  pix  resolution.  Due  to  the  large  imaging  field  required,  60  mm  f/2.8  Nikkor  lenses 
were  used.  To  reduce  reflections,  the  model  was  painted  mat-black  and  monochromatic  filters,  with  a 
corresponding  wavelength  ofX  =  532  nm,  were  installed  on  the  lenses.  The  light  sheet  was  approximately 
2  mm  in  thickness,  set  parallel  to  the  flow  direction  and  aligned  onto  the  airfoil  quarter-span  position.  Due  to 
the  large  imaging  area  of  0.0864  m^  the  laser  power  was  set  to  90%  for  both  cavities.  With  the  use  of 
compressed  air  driven  through  four  Laskin  nozzles,  DEHS  seeding  particles  less  than  1  pm  in  diameter  were 
introduced  into  the  settling  chamber  using  a  vertical  rake  aligned  with  the  measurement  plane. 

The  oscillation  rig  consists  of  a  base  structure,  a  set  of  linear-motors  connected  with  each  other  via  a 
linkage  system  and  a  carbon-fibre  SD7003  wall-spanning  profile  weighing  306  g.  The  profile  has  a  chord 
length  of  120  mm  and  a  span  of  450  mm.  The  profile-tip  spacing  at  the  walls  is  less  than  2  mm  on  either 
side.  Maximum  static  blockage  in  the  test  section  based  on  the  frontal  area  was  under  3%.  A  schematic  of 
the  tandem  configuration  integrated  into  the  Eiffel-type  wind  tunnel  test  section  at  the  Institute  of  Fluid 
Mechanics  and  Aerodynamics  (TU  Darmstadt)  is  shown  in  Figure  2-9.  The  various  components  of  the 
standard  Particle  Image  Velocimetry  (PIV)  system  used  for  these  measurements  are  also  included  in  the 
schematic. 
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Figure  2-9:  Test  Section  of  TU  Darmstadt  Wind  Tunnei,  Showing  Two  Linear  Motor  Pairs  for 
Tandem-Wing  Pitch-Piunge;  Oniy  One  Set  of  Motors  and  One  Wing  was  Used  in  the 
Present  Study:  (a)  PiV  Cameras;  (b)  Beam  Expander;  (c)  PiV  image  Frames; 

(d)  Waii-Spanning  Carbon-Fibre  Profiie;  (e)  Embedded  Piezo-Eiectric 
Force  Sensors;  (f)  Test  Section;  (g)  Laser  Head;  (h)  Linear 
Motors  with  Linkage  System;  and  (i)  Base  Structure. 


LinMot  PSOl-48  x  240F-C  linear  motors  used  to  drive  the  pitch-plunge  motion.  A  displacement  accuracy  of 
<0.5  mm  was  achieved  in  the  setup,  depending  on  the  duration  of  the  experiment  and  the  temperature  of  the 
bushings.  The  control  of  the  linear  motors  was  achieved  with  a  combination  of  digital  inputs  for  high-  speed 
triggering  and  a  serial  port  communication.  A  LabView  8.2  control  program  was  developed  to  operate  the 
communication  system  with  the  linear  motors.  Setup  and  configuration  of  the  motors  was  done  via  LinMot’s 
own  configuration  and  test  program  (LinTalk).  Additional  external  position  sensors  were  mounted  on  the 
motor  units  for  higher  positional  accuracy,  allowing  for  a  dynamic  angle-of-attack  accuracy  of  less  than 
0.5  deg. 

PIV  image  pairs  were  sampled  at  15  Hz  allowing  for  6  phases  to  be  recorded  per  cycle  at  ^  =  0.25. 
In  order  to  construct  the  ensemble  velocity  fields  of  12  phases  per  cycle,  two  staggered  sets  with 
100  images  per  phase  were  ensemble-averaged.  In  all  cases  the  first  two  starting  cycles  were  removed 
from  all  ensembles.  Each  camera  imaged  a  field  corresponding  to  x/c  =  2  and  y/c  =  1.5,  with  a  resolution 
of  800  pix/c  (6.7  pix/mm).  Reflections  on  the  model  surface  were  strongest  at  the  bottom  of  the  stroke 
where  a  region  0.04  c  normal  to  the  airfoil  surface  was  deemed  to  be  unreliable.  Shadows  and  strong 
reflections  on  the  pressure  (lower)  side  required  masking.  Parallax  effects  were  strongest  at  the  top  of  the 
stroke  and  at  this  position  were  responsible  for  hiding  a  region  0.03  c  normal  to  the  airfoil  surface. 
The  vector  fields  were  calculated  using  an  adaptive  correlation  with  32  x  32  pix  interrogation  windows 
and  a  50%  overlap.  A  3  x  3  filter  was  used  to  lightly  smooth  the  vector  fields  in  order  to  more  clearly 
define  the  vortical  structures  in  the  wake.  A  local  neighbourhood  validation  using  a  9  x  9  moving-average 
filter  and  an  acceptance  factor  of  0.2  was  employed  to  eliminate  outliers.  This,  however,  also  had  the 
effect  of  smoothing  the  velocity  gradients,  thereby  thickening  the  shear  layers. 
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2.3.4  University  of  Florida  Wind  Tunnel 

The  UF-REEF/AFRL  low  Reynolds  number  Aerodynamic  Characterization  Facility  (ACF)  is  an  open  jet 
wind  tunnel.  The  test  section  is  a  42”  square  column  of  air  with  a  length  of  15  ft.  What  is  unique  about  the 
ACF  is  that  it  offers  uniform  air  flow  at  mean  free  stream  velocities  ranging  from  1  to  22  m/s.  Turbulence 
intensities  measured  are  less  than  0.11%  between  free  stream  velocities  of  1  and  7.5  m/s.  Photographs  of 
the  tunnel  are  given  in  Figure  2-10. 


Figure  2-10:  University  of  Fiorida  REEF  Open-Jet  Wind  Tunnei;  Exterior  View  (ieft) 
and  Diffuser  Beiimouth  as  Viewed  inside  Chamber  Enciosing  Test  Section  (right). 

The  ACF’s  Dynamic  Pitching-Plunging  Motion  Rig  (DPPMR)  consists  of  two  Parker  ironless  linear  motors 
each  connected  to  an  Aries  model  AR-20AE  driver.  Each  linear  motor  has  a  767  mm  travel  length  with  an 
encoder  of  5  micron  resolution.  The  Aries  drivers  are  controlled  by  a  Galil  DMC-2020  motion  controller. 
Based  on  a  user  input  through  a  LabView  interface,  the  rig  can  be  commanded  to  output  a  standard  TTL 
trigger  at  a  desired  encoder  position,  allowing  for  reliable  PIV  phase  averaging. 

A  LaVision  particle  image  velocimetry  (PIV)  system  was  used  to  make  velocity  field  measurements, 
typically  at  four  instantaneous  phases  (0°,  90°,  180°,  and  270)  throughout  the  cyclical  motion.  The  laser  is  a 
Litron  Nano  L  135-15,  double-pulsing  at  135  mJ/pulse  perpendicular  to  the  model  planform.  The  camera  is 
an  Imager  ProX-4M  (2048  x  2048  pixel  CCD)  capable  of  acquiring  images  at  14  If  ames/sec  for  an  external 
trigger.  A  Nikon  60  mm  lens  with  an  f-number  of  2.3  was  used.  These  components  were  controlled  through  a 
system  computer  running  LaVision  DaVis  software.  The  system  computer  was  externally  triggered  from  the 
TTL  signal  from  the  DPPMR.  This  signal  initialized  the  laser/camera  trigger  sequence.  A  time  delay  of 
300  ps  was  assigned  between  the  image  pairs  of  the  camera  which  resulted  in  19  pixel  displacement  in  the 
free  stream.  The  Imager  ProX-4M  was  camera  oriented  to  view  a  two  dimensional  cross-section  of  the  fiat 
plate  at  75%  of  the  span.  For  all  of  the  measurements  presented,  the  streamwise  domain  of  interest  was  split 
into  two  domains,  each  approximately  119  mm  long.  This  allowed  for  increased  resolution  with  resulting 
grids  1.87  mm  square.  The  flow  was  seeded  with  a  Laskin  nozzle  seeder  using  Olive  Oil  as  the  medium. 
This  seeder  generated  particles  approximately  1  pm  in  diameter  at  a  rate  of  7  x  1010  particles  per  second. 

DaVis  7.2  PIV  software  package  from  LaVision  was  used  to  analyze  the  PIV  images  using  a  multi-pass 
cross-correlation  method.  The  first  pass  consisted  of  an  interrogation  region  of  64  x  64  pixels  with  50% 
overlap.  Two  more  passes  were  conducted  using  a  32  x  32  pixel  interrogation  window  with  50%  overlap. 
As  a  rough  filtering  mechanism,  DaVis  was  used  to  post  process  each  single  vector  field  to  remove  vectors 
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which  are  greater  than  2.5  times  the  root  mean  squared  (rms)  value  of  its  neighbours.  This  procedure  is 
completed  for  420  images  at  a  single  phase  of  motion.  As  a  final  processing  step,  the  mean  velocity  field  is 
calculated  by  averaging  each  interrogation  region  through  all  the  valid  velocity  vectors  throughout  all  the 
ensembles.  Any  velocity  vector  that  resulted  in  a  deviation  outside  of  two  standard  deviations  throughout 
the  ensemble  was  also  removed.  This  effectively  removes  any  erroneous  velocity  vectors  calculated  from 
lack  of  seeding  between  images. 


2.4  FORCE  BALANCES 

Force  measurement  in  unsteady  aerodynamics  is  notoriously  difficult.  Especially  in  wind  tunnels,  where 
physical  motion  rates  are  high,  subtraction  of  inertial  loads  is  a  challenge,  because  inertial  loads  can  easily 
overwhelm  the  aerodynamic  loads  that  one  ultimately  wishes  to  extract.  Sting  and  model  vibration  produce 
additional  force  errors  and  cycle-to-cycle  variations.  Aggressive  low-pass  filtering  is  generally  required, 
but  of  course  one  must  take  care  that  the  filter  does  not  remove  too  much  of  the  “true”  aerodynamic 
component.  These  difficulties  are  quite  attenuated  in  water-based  facilities,  where  physical  motion  rates  are 
much  smaller,  but  there  one  must  contend  with  the  problem  of  fluid  conductivity  and  rusting  or  oxidation  of 
metals.  Thus,  a  large  part  of  the  present  work  was  the  development,  calibration  and  proof-of-concept  testing 
of  force  balances. 

2.4.1  AFRL  Fibre-Bragg  Grating  Balance 

To  circumvent  difficulties  associated  with  electrical  strain  gauges  in  water  -  waterproofing,  routing  of 
wires,  drift,  gauge  delamination  and  so  forth  -  an  optical  approach,  using  fibre-Bragg  gratings  (FBGs) 
[60],  was  used.  Coherent  light  is  sent  through  the  fibre  and  through  gratings  written  onto  the  fibre. 
Each  grating  reflects  light  of  very  narrow  bandwidth.  If  the  segment  of  the  fibre  containing  a  grating  is 
strained,  the  reflected  light  wavelength  shifts  proportionately.  Strain  of  the  fibre  could  be  due  to 
mechanical  strain  of  the  underlying  substrate  (the  flexure  joint  in  the  force  balance)  and  to  thermal  effects, 
which  must  be  removed  though  appropriate  compensation.  In  the  present  application,  a  single  fibre  with 
4  grating  elements  was  integrated  into  a  two-flexure-joint  airfoil  mount,  thus  serving  as  an  integrated  force 
balance  (Figure  2-11).  The  balance  can  resolve  axial  force,  normal  force  and  pitching  moment,  though 
only  the  lift  is  reported  here. 


Figure  2-11:  3-Component  Force  Balance  Based  on  Fibre  Bragg  Grating  (FBG)  Sensors, 
Integrated  with  Airfoil  Mount;  Schematic  (left)  and  Photograph  from 
Underside,  with  Detail  of  a  Flexure  Joint  (right). 


The  FBG  signal  was  interrogated  via  a  Micron  Optics  sml30  instrument  [61],  with  sampling  at  250  Hz 
phase-averaged  over  200  periods  of  oscillation  and  low-pass  filtered  at  6.5  Hz,  with  the  first  5  periods 
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removed  from  the  sample.  Inertial  tares  were  conducted  by  draining  the  water  tunnel  and  repeating  the 
airfoil  motion. 

The  balance  was  calibrated  in  benchtop  tests  with  loads  applied  at  0°  and  ±45°.  These  data  were  used  to 
compute  the  load  cell  calibration  matrix.  The  standard  error  of  the  calibration  matrix  for  the  lift  force  is 
0.16  N,  which  corresponds  to  a  lift  coefficient  standard  error  of  0.03. 

2.4.2  Cambridge  University  External  Balance 

Two-component  force  balance  measurements  were  taken  in  both  air  and  water  for  the  wing  waving 
through  90  degrees  or  5.37  chord-lengths  of  travel  at  3/4  span.  Data  was  sampled  at  7  kHz  with  and 
without  a  30  Hz  low  pass  filter.  Wing  lift  and  drag  coefficients  were  obtained  by  subtracting  the  inertial 
forces  measured  in  air  from  the  forces  measured  in  water  and  normalizing  by  the  local  wing  velocity  such 
that  Cl=  6L/[\rho\omega^2c(r_t^3-r_r^3)]  where  r_t  and  r_r  are  the  distances  from  the  axis  of  rotation  to 
the  wing  tip  and  wing  root.  Force  data  was  averaged  over  five  runs. 

The  force  balance  is  capable  of  measuring  two  force  components  to  10  N  with  a  resolution  of  at  least 
0.01  N.  When  paired  with  a  Fylde  FE-379-TA  transducer  amplifier  and  National  Instruments  USB-622I 
Multi-Function  DAQ,  unsteady  lift  and  drag  forces  can  be  measured  at  frequencies  up  to  250  kHz. 

Figure  2-12  shows  the  lift  coefficient  values  calculated  from  the  unfiltered  data,  data  acquired  using  a 
hardware  30  Hz  low  pass  filter,  and  the  same  data  with  a  moving-average  smoothing  applied. 
The  unfiltered  lift  signal  contains  high-frequency  electrical  noise  near  100  Hz  and  a  mechanical  vibration 
between  15  and  20  Hz.  This  vibration  is  evident  in  both  air  and  water  when  the  motor  is  run  with  and 
without  the  wing  attached.  There  is  frequency  shift  between  air  and  water  thus  the  inertial  data  acquired  in 
air  is  smoothed  before  it  is  subtracted  from  the  force  data  acquired  in  water.  The  force  data  shown  in  the 
following  plots  was  acquired  with  the  low  pass  filter  in  place  to  remove  the  electrical  noise,  the  inertial 
data  subtracted,  then  smoothed  to  eliminate  the  signal  from  the  vibrations. 
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Figure  2-12:  Sampled  CUED  Tow  Tank  Force  Balance  Data,  with  Various  Types  of  Signal  Filtering. 
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2.4.3  TU  Braunschweig  Internal  Balance 

The  3D  model  support  of  TU  BS  is  designed  for  high  precision  movements  with  3  degrees  of  freedom  in 
the  upright  plane.  The  model  support  is  driven  by  linear  synchronous  direct  drives  with  moving  coil, 
with  a  maximum  of  700  N  or  a  speed  of  5.3  m/s,  during  high  dynamic  and  high  precision  movements. 
The  standard  movements  are  a  plunging  motion  up  to  ±150  mm  @  3  Hz  and  a  pitching  motion  with  max. 
±20°  @  3  Hz.  The  position  accuracy  is  ±0.1  mm  or  ±0.1°:  The  repetitive  accuracy  of  both  types  of  motion 
is  ±0.05  mm  or  ±0.1°. 

A  sketch  of  the  actuation  principle  is  shown  in  Figure  2-13.  The  two  vertical  actuators  are  used  to  create 
pitch  and  plunge  motion  of  the  model.  This  set  up  generally  causes  a  small  longitudinal  motion  in  the 
streamwise  direction.  To  remove  this  problem,  additional  horizontal  actuators  will  be  added  in  the  future. 


Figure  2-13:  Motion  Scheme  of  3D  Model  Support  (left)  and  Photo  of  AR  =  2 
Flat-Plate  Model  Mounted  in  the  TUBS  Wind  Tunnel  Test  Section  (right). 


A  small  six-component  internal  strain  gauge  balance  is  used  to  determine  the  aerodynamic  forces. 
Dynamic  measurements  with  wind-off  conditions  are  performed  to  eliminate  inertial  forces.  Wind-off 
measurements  with  both  the  original  flat-plate  model  and  a  cylindrical  dummy  model  with  the  same  mass 
and  inertia  moments  revealed  a  significant  effect  of  wind-off  aerodynamic  forces  on  the  results.  Hence  the 
final  force  data  is  based  on  the  dummy  model  wind-off  data.  The  wind-on  and  wind-off  data  are  then 
subtracted,  after  averaging  the  raw  data  over  80  periods.  After  this  averaging  there  still  appear  high- 
frequency  perturbations  caused  by  the  harmonics  of  the  natural  vibration  of  the  actuator  mechanism.  These 
are  removed  by  low-pass  filtering  using  a  limit  frequency  of  15  Hz.  More  details  are  found  in  Wokoeck 
et  al.  [62]  and  Moeller  et  al.  [63]. 

2.4.4  ONERA  Force  Balance 

An  internal  balance  is  fitted  inside  the  hollow  wing  model.  For  a  10%  thick  airfoil,  the  minimum  chord 
necessary  to  fit  the  balance  is  12.5  cm.  The  centre  of  the  balance  is  placed  15  cm  under  the  free  surface  to 
improve  its  accuracy;  but  the  wing  has  to  pierce  the  free  surface  and  therefore  the  balance  housing  must  be 
watertight.  A  linear  potentiometer  placed  directly  on  the  wing  axis  indicated  that  for  a  desired  pure-plunge 
motion,  there  was  a  pitch  oscillation  of  ~  0.5°. 

As  the  towing  tank  can  be  easily  filled  or  emptied,  tares  without  water  were  run  to  analyze  the  electrical 
signals  of  the  balance.  The  only  significant  forces  generated  during  sinusoidal  translations  in  horizontal 
movements  in  air  are  due  to  the  wing  inertia.  The  weight  of  the  wing  is  well  known  but  the  inertial  forces  are 
a  little  higher  than  the  product  of  its  mass  by  its  acceleration,  because  a  part  of  the  balance  weight  is 
measured.  The  important  variations  of  the  electric  signals  are  due  to  the  vibrations  induced  by  the  oscillation 
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mechanism  and  by  the  towing  tank  carriage  system.  Tests  are  run  with  and  without  motion  of  the  carriage 
system.  Rubber  isolators  are  placed  between  the  oscillation  system  and  the  frame  to  avoid  transmission  of 
the  carriage  vibrations.  The  level  of  vibration  with  and  without  the  longitudinal  velocity  is  the  same; 
the  structural  frequency  of  the  frame  is  >10  Hz. 

It  appears  that  the  structural  frequency  of  the  wing  is  of  about  8  Hz  in  the  air  (Figure  2-14).  Data  is  acquired 
at  2  kHz  without  any  filter.  Numerical  filtering  is  adjusted  to  maintain,  for  every  periodic  plunging  motion, 
some  small  fluctuations  of  the  signal.  The  value  of  the  mass  to  be  considered,  such  that  the  RMS  value  of  the 
measured  minus  the  inertial  forces  is  minima,  is  computed  from  several  tests  in  air.  The  residual  forces  are 
less  than  0.1  N;  that  is  small  compared  to  the  forces  of  about  10  N  measured  in  the  water.  In  the  water  the 
vibration  of  the  wing  is  rapidly  damped  and  the  structural  frequency  is  decreased  (3.5  Hz)  due  to  the  virtual 
mass  effect.  A  compromise  has  to  be  found  between  the  balance  sensitivity  and  the  rigidity  necessary  to  get  a 
high  structural  frequency. 


static  AIR 

i=8°  SD7003  AR=14  V=0  h=0.5  f=0.51  Hz 


Figure  2-14:  Inertial  Tares  in  Air,  ONERA  Tow  Tank. 


2.4.5  TU  Darmstadt  Force  Balance 

For  direct  lift  and  pitching  moment  measurements  in  the  TU  Darmstadt  wind  tunnel,  a  pair  of  one- 
component  Kistler  9217A  piezo-electric  force  sensors  were  integrated  directly  below  the  SD7003  profile, 
one  at  the  quarter-chord  position  (point  of  rotation)  and  one  at  the  trailing  edge,  as  shown  in  Figure  2-15. 
Together  the  two  sensors  measure  the  inertial  and  aerodynamic  forces  during  the  prescribed  movement. 
Similar  piezo-electric  sensors  have  been  used  successfully  for  the  measurement  of  pitching  and  plunging 
foils  in  a  tow-tank,  both  for  a  single  foil  and  in  the  wake  of  a  cylinder;  see  Anderson  et  al.  [18]  and  Beal 
et  al.  [64],  respectively.  The  total  static  tare  weight  of  the  system  (airfoil  and  linkage  combined)  was 
394  g.  The  analogue-output  charge  signals  from  the  two  piezo-electric  force  sensors  were  sent  through  the 
wind-tunnel  floor  to  a  Kistler  5073A411  charge  amplifier,  which  in  turn  converted  the  signals  into  an 
analogue  voltage  which  was  fed  into  a  16-bit  National  Instruments  6259  A/D  board.  Since  the  physical 
and  sting  natural  frequencies  were  on  the  order  of  1  Hz  and  20  Hz,  respectively,  the  signals  were  run 
through  a  5  Hz  low-pass  filter  in  LabView  8.2  and  then  further  post-processed  using  MATLAB  7.3. 
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Figure  2-15:  Location  of  Piezo-Eiectric  Force  Sensors  Directiy  Beiow  SD7003  Profiie; 

Note  Sensor  1  Located  at  Quarter  Chord  Position  (Point  of  Rotation) 
and  Sensor  2  Located  at  Trailing  Edge. 

All  lift  and  moment  measurements  were  based  on  an  ensemble  of  30  clean  cycles,  sampled  at  1  kHz, 
where  the  first  four  cycles  as  well  as  the  last  cycle  were  discarded  to  avoid  aerodynamic  and  inertial 
starting  and  stopping  effects.  To  subtract  the  inertial  tare,  a  corresponding  ensemble  of  30  clean  cycles 
was  measured  with  the  wind  tunnel  turned  off  This  technique  proved  to  be  very  repeatable  (<  1%)  for  the 
lower  reduced  frequencies  (k  <0.15)  where  the  aerodynamic  contribution  was  always  equal  or  larger  than 
the  inertial  contribution.  Results  at  these  lower  reduced  frequencies  are  presented  in  Rival  and  Tropea 
[65],  where  the  accuracy  was  estimated  to  lie  within  ACi  =  ±0.05  and  ACm=  ±0.02,  respectively.  However, 
at  higher  reduced  frequencies  (k  >  0.2),  the  inertial  loads  grew  with  the  frequency  squared  and  therefore  very 
quickly  dominated  the  total  measured  forces,  thus  degrading  the  accuracy  of  the  measurements.  For  these 
cases  a  dummy  cylinder  was  used  to  determine  the  wind-off  aerodynamic  contribution  (Figure  2-16)  but  had 
only  limited  success.  For  more  details,  see  Rival  et  al.  [66]. 


Figure  2-16:  Positioning  of  Dummy  Cylinder  at  Center-of-Mass  Position  Inside  Wind-Tunnel  Test 
Section  (left)  and  Comparison  of  Dummy  Cylinder  with  Original  SD7003  Profile; 

Note  Cylinder  Identical  to  Profile  in  Mass  (±0.5  g)  and  Spanwise  Distribution. 
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2.5  COMPUTATIONAL  METHODS 

Computations  range  from  analytical  formulas  based  on  conformal  mapping  and  inviscid  2D  vortex  particle 
methods,  to  commercial  and  custom  in-house  RANS  methods,  to  higher-order  3D  LES  methods.  Each  has  its 
place  and  each  has  tradeoffs  of  computational  cost/complexity  of  implementation  vs.  accuracy.  One  of  our 
objectives  is  to  assess  what  level  of  accuracy  is  necessary  for  this  or  that  objective,  such  as  lift  coefficient 
prediction  (at  a  minimum)  to  capturing  the  main  flowfield  features  such  as  shed  vortices,  all  the  way  to 
accurately  resolving  boundary  layer  transition. 


2.5.1  AFRL  Implicit  Large  Eddy  Simulation 

The  governing  equations  are  the  unfiltered  full  compressible  Navier-Stokes  equations  cast  in  strong 
conservative  form  after  introducing  a  general  time-dependent  curvilinear  coordinate  transformation 
(x,  y,  X,  t)  —>  rj,  Q  t)  from  physical  to  computational  space.  In  terms  of  non-dimensional  variables, 
these  equations  can  be  written  in  vector  notation  as: 


dt 


+ 


OH, 


J_ 

Re 


drj  d(^ 


the  vector  of  dependent  variables  and  inviscid  fluxes  are  given  by: 


pu  pv  pw 


pE,]' 


(1) 


(2) 


pU 

pV 

pW 

puU  +  ^^p 

puV  +  ri^p 

•"'=7 

P^W  +  C^p 

pvU  +  ^^p 

pvV  +  r]yp 

pvW  +  ^^p 

pwU  +  ^^p 

pwV  +  ri^p 

pwW  +  ^^p 

_pEp  +  pU_ 

pE,V  +  pV 

pE,W  +  pW 

with: 


V  =  ri^+ri^u  +  riyV  +  ri^w  =  ri^+V 

E.  —  - r — —  H —  In 


+  -\U^  +V^ 


(4) 


Here,  S^/Sx  with  similar  definitions  for  the  other  metric  quantities.  The  viscous  fluxes,  F^,  G  v 

and  Hy  can  be  found,  for  instance,  in  Anderson  et  al.  [67].  In  the  expressions  above,  u,  v,  w  are  the 
Cartesian  velocity  components,  p  the  density,  p  the  pressure,  and  T  the  temperature.  The  perfect  gas 
relationship  p  =  pT/  (yM^),  as  well  as  Sutherland’s  law  for  viscosity  are  also  assumed.  All  flow  variables 
have  been  normalized  by  their  respective  reference  freestream  values  except  for  pressure  which  has  been 
normalized  by  . 
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The  above  governing  equations  correspond  to  the  original  unfiltered  Navier-Stokes  equations,  and  are 
used  without  change  in  laminar,  transitional  or  fully  turbulent  regions  of  the  flow.  Unlike  the  standard 
LES  approach,  no  additional  sub-grid  stress  (SGS)  and  heat  flux  terms  are  appended.  Instead,  a  high-order 
low-pass  filter  operator  (to  be  described  later)  is  applied  to  the  conserved  dependent  variables  during  the 
solution  of  the  standard  Navier-Stokes  equations.  This  highly-discriminating  filter  selectively  damps  only 
the  evolving  poorly-resolved  high-frequency  content  of  the  solution.  This  filtering  regularization 
procedure  provides  an  attractive  alternative  to  the  use  of  standard  SGS  models,  and  has  been  found  to 
yield  suitable  results  for  several  canonical  turbulent  flows  (Visbal  et  al.  [68]  and  Visbal  et  al.  [69])  on 
LES-level  grids. 

All  simulations  are  performed  with  the  extensively  validated  high-order  FDL3DI  Navier-Stokes  solver, 
described  in  more  detail  in  Visbal  and  Gaitonde  [70].  In  this  code,  a  finite-difference  approach  is  employed 
to  discretize  the  governing  equations,  and  all  spatial  derivatives  are  obtained  employing  a  6*-order  compact- 
differencing  scheme.  In  order  to  eliminate  high-frequency  spurious  components,  an  8^^-order  Pade-type  low- 
pass  spatial  filtering  operator  [70]  is  applied  to  the  conserved  variables  along  each  transformed  coordinate 
direction  after  each  time  step  or  sub-iteration.  For  transitional  and  turbulent  flows,  this  filtering  technique 
provides  an  effective  implicit  LES  approach. 

For  the  case  of  a  manoeuvring  airfoil,  the  grid  is  moved  in  rigidly  using  the  prescribed  airfoil  motion. 
To  ensure  that  the  Geometric  Conservation  Law  (GCL)  is  satisfied,  the  time  metric  terms  are  evaluated  as 
in  Visbal  and  Gaitonde  [71].  The  original  airfoil  sharp  trailing  edge  was  rounded  with  a  circular  arc  of 
radius  r/c-  0.0004  to  use  an  O-mesh  topology.  The  computational  mesh  consisted  of  (649  x  395  x  101) 
points  in  the  streamwise,  normal  and  spanwise  direction,  respectively.  Grid  points  were  concentrated  near 
the  airfoil  to  capture  the  transition  process.  For  the  three-dimensional  simulations,  which  invoked 
periodicity  in  the  spanwise  direction,  the  mesh  had  a  span  s/c  =  0.2.  Further  details  on  the  mesh  and 
boundary  conditions  can  be  found  in  Visbal  [98]. 

Plunging  simulations  were  started  from  previously  computed  static  solutions  at  the  corresponding  mean 
angle  of  attack.  Simulations  were  then  advanced  in  time  for  more  than  25  cycles  in  order  to  guarantee  a 
time-asymptotic  nearly-periodic  state.  A  very  small  computational  nominal  time  step  At  U  /c  =  0.00005 
was  prescribed  in  order  to  provide  sufficient  temporal  resolution  of  the  abrupt  spanwise  breakdown  of  the 
leading-edge  vortices.  This  value  of  Zl  ^  corresponds  to  16,000  time  steps  per  cycle  for  a  reduced  frequency 
k=  3.93.  Finally,  all  computations  were  performed  employing  a  low  freestream  Mach  number  M  =  0.1, 
as  required  with  the  present  compressible  Navier-Stokes  solver.  The  good  agreement  of  the  computed  lift 
coefficient  with  the  inviscid  incompressible  theory  for  ^  =  3.93  demonstrated  that  compressibility  effects 
were  minor. 

2.5.2  AFRL  Vortex-Particle  Method 

The  method  presented  in  this  work  is  based  on  a  2D  inviscid  point  vortex  method  similar  to  those 
presented  in  Katz  [72]  and  Karamcheti  [73].  Several  important  modifications  have  been  incorporated. 
First,  conservation  of  circulation  is  implicitly  satisfied,  rather  than  being  imposed  as  an  additional  constraint. 
Second,  a  vortex  decay  model  is  made  available  for  flow  that  is  separated  into  the  wake. 

In  two  dimensions,  a  body  can  be  represented  by  a  collection  of  bound  vortices  particles.  Each  vortex 
particle  can  be  interpreted  as  an  infinitely-long  vortex  filament  oriented  normal  to  the  plane.  Each  vortex 
induces  a  purely  rotational  velocity  field.  According  to  the  Biot-Savart  law,  the  magnitude  of  the  induced 
velocity  varies  directly  with  the  strength  of  the  vortex,  F,  and  inversely  with  the  distance  between  the 
vortex  and  the  point  of  interest.  Since  this  results  in  an  infinite  velocity  at  the  location  of  a  vortex  particle, 
in  practice  the  velocity  must  be  kept  finite  according  to  some  rule.  Although  several  possibilities  exist, 
we  use  the  following  formulation  based  on  the  error  function: 
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=  —  - - jerf  I  * -  I  r  X 

2.T|r-|=‘ 


(5) 


where  v  is  the  velocity  vector,  f  is  displacement  vector  from  the  vortex  to  the  point  of  interest,  k  is  the 
unit  normal  vector  out  of  the  plane,  and  r cutoff the  cutoff  radius  below  which  the  inverse  relationship  with 
r  does  not  hold. 


For  a  thin  body  represented  by  Nbv  bound  vortices  ordered  along  the  body’s  length,  we  solve  for  the 
strengths  of  the  bound  vortices  at  a  given  moment  in  time.  We  introduce  a  control  point  in  between  each 
pair  of  bound  vortices.  Let  Ncp  =  Nbv  “  1  be  the  number  of  control  points.  At  each  control  point,  assign  a 
strength  representing  the  strength  of  a  ring  of  vorticity.  In  2D,  we  imagine  the  ring  to  extend  to  ±oo  in  the 
out-of-plane  direction,  so  the  two  segments  parallel  to  the  2D  plane  have  no  contribution  to  the  flow. 

The  strength  of  the  bound  vortices  T  are  related  to  the  control  point  strengths  /  by: 


f 


=  T1  +  S, 


(6) 


where  T  is  a  mapping  matrix  with  elements  of +1  on  the  main  diagonal  and  -1  on  the  lower  off-diagonal. 

S  is  the  amount  of  vorticity  that  has  been  shed  into  the  wake  from  the  given  bound  vortex  at  the  previous 
moment  in  time,  if  any.  In  this  manner,  vorticity  is  implicitly  conserved  since  each  vortex  ring  contributes 
a  net  vorticity  of  zero.  At  each  time  step,  the  velocity  at  each  control  point  is  calculated  by  enforcing  flow 
tangency.  Let  /  be  an  Ncp  x  Nbv  influence  matrix  with  elements  iij  equal  to  the  velocity  induced  by  a 

bound  vortex  j  of  unit  strength  on  control  point  /.  Let  be  a  vector  of  length  Ncp  with  elements  v .  that 

are  the  normal  component  of  velocity  at  the  bound  vortex  /,  equal  to  the  velocity  of  the  body  plus  the 
velocity  induced  by  the  bound  vorticity  and  the  wake.  Flow  tangency  leads  to: 


/f=  -.V 

ITl  =  -IS  -  N  (7) 

Once  the  bound  vortex  strengths  have  been  computed  for  a  given  time  step,  vorticity  is  shed  into  the  wake 
at  the  trailing  edge.  A  new  wake  panel  is  created  a  distance  behind  the  trailing  edge  panel  equal  to  the  time 

step  times  the  local  velocity.  The  strength  of  this  panel  is  set  equal  to  i?  where  k  is  an 

attenuation  factor  typically  taken  to  be  one.  The  bound  vortex  strength  and  the  wake  vortex 

strength  can  then  be  calculated.  Once  the  vortices  have  been  shed  and  convected,  the  state  is  ready 

to  be  advanced  by  another  time  step. 

The  resulting  calculation  gives  values  of  lift  coefficient  time  history,  accounting  for  airfoil  camber  (but  not 
explicitly  for  thickness)  and  time-history  of  vorticity  shed  form  the  trailing  edge.  It  does  not  however 
account  for  leading-edge  vortex  formation  or  shedding. 


2.5.3  METU  RANS  Approach 

Computations  use  Fluent  [74]  -  an  unsteady,  pressure  based  Reynolds-Averaged  Navier-Stokes  solver 
with  second  order  upwind  spatial  discretization.  The  pressure-velocity  coupling  is  handled  by  the  SIMPLE 
algorithm.  Motion  kinematics  is  implemented  using  ‘dynamic  mesh’  feature  of  Fluent.  To  model  the 
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motion,  the  whole  grid  is  divided  into  two  parts:  inner  grid  with  20  c  diameter  and  outer  grid  with  60  c 
diameter.  In  solution  process,  the  inner  grid  around  the  airfoil  performs  the  prescribed  motion  while  the 
outer  grid  is  stationary  and  deforming  with  appropriate  rules.  The  dynamic  mesh  feature  of  the  code  limits 
the  unsteady  formulation  to  the  first  order  in  time  and  so  for  all  calculations  first-order  temporal 
discretization  is  used.  The  far-field  boundary  conditions  are  set  as  velocity  inlet  and  pressure  outlet. 
The  free-stream  flow  is  assumed  to  be  turbulent  and  Menter’s  Shear  Stress  Transport  model  is  used  for 
turbulence  closure,  which  is  a  favourable  model  in  prediction  of  adverse  pressure  gradients  [75]. 

An  unstructured  O-type  grid  with  outer  diameter  of  60  c  is  used.  The  y+  value  for  all  cases  is  chosen  as 
less  than  1.  Figure  2-17  shows  the  whole  domain  and  the  grid  around  the  airfoil.  A  grid  refinement  process 
is  conducted  to  have  grid-independent  solution.  The  previous  studies  showed  that  a  grid  with  92358  cells 
is  enough  to  have  the  grid-independent  solution  for  SD  7003  airfoil  [76].  For  flat  plate  with  2.5%  thickness, 
132804  cells  are  used.  Moreover,  960  time  steps  per  period  are  chosen  for  all  calculations.  The  results  are 
represented  for  [6T-7T]  period. 


Figure  2-17:  Computational  Mesh  for  SD7003  Airfoil; 
Full  Grid  (left)  and  Zoomed-ln  to  Airfoil  Region  (right). 


2.5.4  University  of  Michigan  RANS  Approach 

Computations  by  the  University  of  Michigan  are  unsteady  RANS,  closed  with  Menter’s  SST  turbulence 
model  [77,78].  The  momentum  equation  becomes: 


the  eddy  viscosity,  Vf,  given  by  5  =  is  the  invariant  measure  of  the  strain  rate. 

Uf  is  the  velocity  component  in  the  direction,  is  the  component  of  the  position  vector,  t  is  time,  p  is 

density,  p  is  pressure,  v  is  the  kinematic  viscosity.  The  transport  equation  for  turbulent  kinetic  energy  k  is: 


dk  d  ^  5  f  ^  dk  ) 


dxj 


The  transport  equation  for  co  is: 


j+  2(1 


1 


dk 

dx^  dx^ 
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Finally,  the  quantities  /?,  /?*,  r,  ^2  ^^e  defined  as  in  Menter’s  SST  formulation  [77,78]. 

Compared  to  Menter’s  original  SST  turbulence  model,  a  limiter  was  built  in  to  the  production  term, 

du[fdui  diEj\ 

in  the  turbulence  kinetic  energy  (TKE)  equation,  as  — v 

&xi\dxj  dxj 

=  minCFjj,  10  ■  where  P^^  is  the  production  term  in  the  original  SST  formulation,  to  prevent  the 

build-up  of  turbulence  in  stagnation  regions.  Another  change  is  the  use  of  invariant  measure  of  the  strain- 
rate  tensor  in  the  formulation  for  the  eddy  viscosity  instead  of  the  vorticity  magnitude,  O  =  ^."’20gQg-. 
The  strain-rate  invariant  is  considered  to  be  a  better  measure  for  the  fluid  deformation,  since  the 
Boussinesq  approximation  is  also  based  on  the  strain-rate.  The  two  differences  between  the  original  and 
the  modified  SST  formulation  are  summarized  in  Table  2-2. 


Table  2-2:  Original  and  Modified  SST  Turbulence  Model,  U  of  Michigan  Computation 


Original  SST 

Modified  SST 

Production  T erm  of  TKE 

duif&Ui  duf\ 

Fjj  =  rainCFjj.  10  - 

Equation 

Eddy  Viscosity 

These  equations  are  solved  with  the  in-house  solver  Loci-STREAM  [79],  a  parallelized  unstructured 
curvilinear  pressure-based  finite-volume  code  with  moving  grid  capabilities.  The  present  calculations  use 
implicit  first  order  time  stepping.  The  convection  terms  are  treated  using  the  second  order  upwind 
scheme  [80,81]  while  pressure  and  viscous  terms  are  treated  using  second  order  schemes.  The  geometric 
conservation  law  [82,83],  a  necessary  consideration  in  domains  with  moving  boundaries,  is  satisfied. 
The  numerical  solutions  were  computed  in  an  open  bounded  domain  with  Loci-STREAM  on  unstructured 
grids  of  different  resolutions.  More  details  can  be  found  in  Kang  et  al.  [84]  and  Kang  et  al.  [85]. 

2.5.5  TU  Darmstadt  RANS  Approach 

Simulations  used  Fluent  6.3.26,  and  were  run  with  second-order  upwind  spatial  and  first-order  temporal 
discretization,  where  the  latter  was  limited  by  the  dynamic-meshing  process  in  the  code.  To  compensate 
for  this  limitation  in  the  temporal  discretization,  fine  time-stepping  of  0.5  ms  was  used,  equivalent  to 
800  time  steps  per  cycle  for  Re  =  30000  and  k  =  0.25.  Further  refinement  of  the  time-stepping  was 
performed  to  check  for  a  time-step  independent  solution  analogous  to  the  grid-independency  tests 
described  below.  20  to  30  iterations  were  run  for  a  given  time  step  with  a  convergence  criterion  based  on 
residuals  dropping  0(10"^)  in  magnitude.  Pressure-velocity  coupling  was  performed  using  the  SIMPLE 
scheme.  For  all  simulations  the  procedure  was  to  run  two  complete  cycles  to  establish  the  appropriate 
initial  conditions  and  then  to  begin  recording  data  on  the  following  cycle. 

The  numerical  domain  was  modelled  after  the  TU  Darmstadt  Eiffel-type  open-return,  closed  test  section 
wind  tunnel.  The  present  computations  were  a  sub-set  of  a  study  on  tandem  airfoils,  but  here  we  report 
only  on  the  single  airfoil;  Figure  2-18  depicts  tandem  airfoils.  The  domain  extends  the  entire  length  of  the 
wind-tunnel  test  section,  and  is  filled  with  block-structured  cells  in  sections  a)  and  d)  before  and  after  the 
dynamic  meshing  region,  as  shown  in  Figure  2-18.  In  these  block-structured  regions  the  skewness  is  zero, 
and  the  aspect  ratio  and  cell  growth  are  nearly  one. 
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Figure  2-18:  Top:  Numerical  Domain  Representative  of  Wind-Tunnel  Test  Section;  Zones  a)  and 
d)  Contain  Pure  Block-Structured  Cells  Whereas  Zones  b)  and  c)  Contain  Dynamic  Layering  and 
Sliding  Zones;  Bottom:  Middle  of  Domain;  Zones  b)  and  c),  with  Two  Circular-Structured  Cores 
Encompassing  the  Airfoils.  Tandem  airfoils  depicted,  but  only  one  airfoil  was  computed  here. 


The  circular-core  zones  surrounding  the  airfoils  provide  for  the  pitching  movement  (through  non-conformal 
boundaries)  whereas  the  plunging  movement  was  realized  through  the  dynamic-layering  method  where  rows 
of  structured  cells  above  and  below  the  circular  cores  were  created  and  destroyed.  The  two  unstructured 
regions  in  the  entire  mesh  were  filled  with  quadrilateral  cells,  providing  the  coupling  of  the  non-conformal 
boundaries  between  the  circular  cores  and  the  block-structured  outer  region.  This  dynamic-meshing 
arrangement  allows  for  high  pitch-  and  plunge  amplitudes  while  maintaining  good  mesh  quality. 

Both  airfoil  and  wind-tunnel  wall  boundary  layers  were  fully-resolved  satisfying  the  y+  =  1  condition 
throughout.  Cell  growth  in  all  boundary  layers  was  limited  to  1.25.  The  trailing  edges  of  the  airfoils  were  set 
with  radii  of  0.1  mm  so  as  to  be  more  representative  of  actual  experimental  airfoils  and  to  ease  the  meshing 
constraints.  An  overall  maximum  skewness  of  0.41  and  maximum  aspect  ratio  of  4.1  existed  within  the 
circular  cores.  A  grid-independency  study  was  performed  and  a  resulting  grid  of  approximately  32000  cells 
was  established.  The  relatively  low  number  of  cells  and  thus  fast  simulation  times  (on  the  order  of  six  hours 
on  a  single  Linux  machine)  can  be  attributed  to  the  nearly  complete  block-structured  mesh. 
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2.5.6  NRC  Large  Eddy  Simulation 

Computations  at  the  Canadian  National  Research  Council  Institute  for  Aerospace  Research  used  the  in-house 
3D  unsteady  incompressible  LES/RANS  code  INSflow  [86].  Recent  numerical  investigations  of  low- 
Reynolds-number  and  flapping-wing  aerodynamics  can  be  found  in  Yuan  et  al.  [87,88,89].  INSflow  uses  the 
integral  form  of  the  conservation  law  for  mass  and  momentum.  A  fully  implicit  second-order  temporal 
differencing  scheme  makes  the  algorithm  stable  for  large  time  steps.  The  discretization  of  the  convective  and 
diffusive  fluxes  was  carried  out  in  a  co-located  variable  arrangement  using  a  finite-volume  approach  that  was 
second-order  accurate  in  space.  The  coupling  of  the  pressure  and  the  velocity  was  handled  using  the 
SIMPLE  algorithm  [90].  The  continuity  equation  was  transformed  into  a  pressure-correction  equation  that 
had  the  same  general  form  as  the  discretized  momentum  equations.  The  use  of  the  co-located  variable 
arrangement  on  non-orthogonal  grids  required  that  the  SIMPLE  algorithm  be  slightly  modified  to  dampen 
numerical  oscillations.  A  pressure-velocity  coupling  method  for  complex  geometries  used  by  Ferziger  and 
Peric  [91]  was  implemented,  where  an  additional  pressure  gradient  term  was  subtracted  from  the  velocity 
value  at  the  surface  of  the  control  volume  to  prevent  non-physical  oscillations.  A  number  of  two-equation 
turbulence  models  and  two  sub-grid-scale  (SGS)  models  were  implemented  for  RANS  and  EES. 

The  calculations  were  performed  on  a  moving  grid  whose  velocity  was  included  in  the  governing  equations 
[86,92]  in  an  inertial  frame  of  reference.  To  avoid  artificial  mass  sources  generated  by  the  grid  velocity, 
a  space  conservation  law  was  introduced  to  ensure  a  fully  conservative  property  in  the  computations, 
as  applied  by  Demirdzic  and  Peric  [93]. 

Two-dimensional  (2D)  simulations  were  performed  for  both  pure  plunging  and  combined  pitching-plunging 
cases.  The  2D  calculations  were  LES-like  since  the  eddy  viscosity  was  added  to  the  diffusion  term. 
Calculations  were  initiated  from  stationary  flow  conditions.  A  C-mesh  with  737  x  65  grid  points  was  used 
for  the  pure  plunging  SD7003  airfoil  case,  an  O-mesh  with  961  x  65  grid  nodes  was  used  for  the  combined 
pitching-plunging  SD7003  airfoil  case,  and  an  O-mesh  with  461  x  65  grid  nodes  was  used  for  the  combined 
pitching-plunging  fiat-plate  case.  The  C-mesh  was  generated  using  the  hyperbolic  method  proposed  by  Barth 
et  al.  [94],  and  the  O-meshes  were  generated  using  an  algebraic  method  similar  to  the  one  of  Peric  [91]. 
In  these  calculations,  the  far  field  was  set  at  10  chords  away  from  the  airfoil  for  the  O-meshes  and  25  chords 
away  from  the  airfoil  the  C-mesh.  The  far  field  was  thus  correctly  assumed  to  be  sufficiently  distanced  not  to 
warrant  any  circulation  corrections.  These  C-  and  O-grids  were  generated  based  on  experience  gained  from 
previous  work  [87,88]  and  ensured  that  the  streamwise  spacing  in  wall  units  satisfied  Ax^  <  100.  In  the 
calculations  of  the  pitching  cases,  the  inner  part  of  the  O-mesh  pitched  in  unison  with  the  airfoil  with  the 
outer  part  fixed.  The  timestep  was  set  at  At  =  771920,  where  T  was  the  period  of  one  pitching  cycle. 
Computations  were  run  for  30  cycles  for  the  2D  plunging  SD7003,  20  cycles  for  the  2D  pitching-plunging 
SD7003,  and  6  cycles  for  the  2D  pitching-plunging  fiat-plate. 

2.5.7  U  of  Toronto  Lumped-Vortex  Model 

An  unsteady,  lumped- vortex  method  is  used  to  compute  2D  and  3D  lift  and  drag.  The  method  is  an 
extension  of  DeLaurier’s  [95]  unsteady  strip-theory  model  that  was  originally  developed  for  the  design  of 
large  flapping- wing  aircraft  in  the  range  of  1.5  m  to  12  m  span. 

The  method  models  all  bound  circulation,  including  any  bound  leading-edge  vortices,  as  being  lumped  at 
the  quarter  chord  with  the  no  flow-through  boundary  condition  imposed  at  the  three-quarter  chord. 
The  method  accounts  for  the  unsteady  wake  by  computing  the  influence  of  wake  vortices  that  are  shed  one 
quarter  chord  aft  of  the  trailing  edge  and  convected  downstream  at  each  time  step.  Apparent  mass 
(non-circulatory  force)  is  computed  separately  and  added  to  the  normal  force.  Viscous  forces  are  modelled 
by  correlating  the  bound  circulation  at  each  spanwise  station  with  a  parameterized  fit  of  the  static  lift-drag 
polar  of  the  airfoil.  Dynamic  stall  is  modelled  by  a  reduction  in  the  bound  circulation  as  the  leading-edge 
vortex  is  shed  off  the  trailing  edge.  2D  results  are  computed  simply  by  using  an  infinite  wing  with  only 
one  spanwise  panel,  as  a  special  case  of  3D  panelling  (Figure  2-19). 
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Figure  2-19  UTIAS  Lumped  Vortex  Method:  Vortex  Panelling  for  a  3D  Rectangular 
Wing  (left)  and  Schematic  of  Lag  Parameters  for  Dynamic  Stall  Model  (right). 

When  performing  a  parametric  fit  on  the  lift-drag  polar  of  the  desired  airfoil,  the  attached  and  stalled  portions 
of  the  data  are  fit  separately.  This  allows  for  predictions  of  lift  and  drag  on  the  airfoil  at  angles  above  the 
static  stall  angle  when  the  airfoil  is  undergoing  dynamic  stall  delay.  The  delay  is  thus  modelled  as  an 
effective  shift  in  the  stall  angle,  Aus,  that  depends  on  the  rate  of  change  of  angle  of  attack  and  the  time  that 
has  passed  since  the  angle  of  attack  exceeded  the  static  stall  angle  (schematic  in  Figure  2-19).  This  is  an 
extension  of  Prouty’s  dynamic  stall  model  for  helicopter  blades  [96]  that  incorporates  the  elements  of  the 
formation  time  for  a  leading-edge  vortex.  The  model  is  currently  employed  only  for  the  2D  cases  and  needs 
to  be  extended  for  the  generalized  3D  case. 
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Chapter  3  -  PURE-PLUNGE  OF  THE  SD7003  AIRFOIL 

3.1  RESUME  OF  EXPERIMENTAL  AND  COMPUTATIONAL  RESULTS  FOR 
THE  FLOWFIELD:  VELOCITY  AND  VORTICITY 

We  begin  by  comparing  the  various  experimental  and  computational  results  for  the  flowfield  about  the 
airfoil.  Experiments  are  particle  image  velocimetry,  and  in  general  are  limited  to  the  suction-side  of  the 
airfoil,  as  the  pressure-side  is  in  the  PIV  light  sheet  shadow.  Computations  of  course  do  not  have  this 
limitation.  For  velocity,  we  use  the  normalized  streamwise  component,  u/U,  as  the  metric  of  choice. 
Vorticity  is  limited  to  the  out-of-plane  component,  and  normalization  is  by  free-stream  velocity,  U, 
and  airfoil  chord,  c.  For  velocity  the  contour  levels  are  0  to  1.5,  while  for  vorticity  they  are  -36  to  +36, 
unless  otherwise  noted  in  the  respective  plot. 

3.1.1  Computations 

2D  and  3D  FES  computations  from  AFRF  are  given  in  Figure  3-1  and  Figure  3-2,  respectively. 
METU  PANS  computations  are  in  Figure  3-3,  NRC  2D  RANS  are  in  Figure  3-4,  and  the  UM  RANS  are  in 
Figure  3-5. 


RTO-TR-AVT-149 


3-1 


PURE-PLUNGE  OF  THE  SD7003  AIRFOIL 


ORGANIZATION 


Figure  3-1:  2D  Computations,  FDL3DI;  Phases  Phi  =  0,  45,  90,  135,  180,  225,  270  and  315. 
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Figure  3-2:  3D  Computations,  FDL3Di;  Phases  Phi  =  0,  45,  90,  135,  180,  225,  270  and  315; 
Phase-Averaged  Veiocity  (ieft  coiumn),  Phase-Averaged  Vorticity 
(middie  coiumn)  and  instantaneous  Vorticity  (right  coiumn). 
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Figure  3-3:  METU  Computations,  Fiuent,  2D;  Phases  Phi  =  0,  90, 120, 150, 180,  210  and  270. 
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Figure  3-4:  lAR  NRC  Computations,  LES;  Phases  Phi  =  0,  90, 180  and  270. 
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Figure  3-5:  University  of  Michigan  Computations;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 


3.1.2  Experiments 

Velocity  and  vorticity  contour  plots  from  PIV  are  presented  as  follows:  AFRL  water  tunnel  Series  #1  is  in 
Figure  3-6  and  Series  #2  is  in  Figure  3-7;  and  the  UM  water  tunnel  results  are  in  Figure  3-8.  We  note  that 
because  the  AFRL  water  tunnel  model  is  mounted  horizontally  across  the  test  section  with  rig  linkages  on 
the  pressure-side  along  the  midspan,  the  natural  PIV  light  sheet  location  is  the  %  span.  But  the  UM  water 
tunnel  installation  is  with  the  wing  hanging  vertically  with  no  rig  linkages  submerged  in  the  water,  and  the 
most  natural  PIV  light  sheet  location  is  the  midspan.  For  the  first  AFRL  campaign,  the  PIV  light  sheet  was 
somewhat  inboard  of  the  %  span,  and  for  the  second  campaign  it  was  somewhat  outboard  of  this  location. 
We  speculate,  barring  further  information  such  as  from  a  desirable  but  impractical  full  3D  computation 
that  includes  the  tunnel  test  section  sidewalls,  bottom  and  free-surface,  that  the  differences  between  the 
two  AFRL  runs  and  between  the  AFRL  runs  and  UM  run  is  3D  variation,  or  spanwise  variation,  which  is 
especially  strong  at  the  bottom  of  the  downstroke. 


3-6 


RTO-TR-AVT-149 


PURE-PLUNGE  OF  THE  SD7003  AIRFOIL 


Figure  3-6:  AFRL  Water  Tunnel,  Entry  #1 ;  Phases  Phi  =  0,  90, 1 80  and  270. 
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Figure  3-8:  University  of  Michigan  PiV;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 


3.2  REYNOLDS  NUMBER  EFFECTS 

Because  of  the  variation  from  facility  to  facility  and  computation  to  computation,  Reynolds  number 
variations  are  difficult  to  verify  in  detail,  as  for  example  a  facility  with  larger  turbulence  intensity  or  more 
vibration  of  the  motion  rig  may  evince  similar  flowfield  features  at  Reynolds  number  X,  to  those  of  a 
quieter  facility  at  Reynolds  number  2X.  One  can,  of  course,  do  a  Reynolds  number  sweep  in  one  facility, 
but  even  there  the  conclusion  is  not  definite  because  of  issues  with  repeatability  from  run  to  run,  which 
may  be  no  smaller  than  the  Reynolds  number  variation  that  one  seeks  to  identify.  Perhaps  a  statistical 
approach  is  required,  where  runs  at  different  Reynolds  numbers  are  repeated  some  number  of  times, 
in  random  order.  This  is  perhaps  possible  for  dye  injection,  but  not  for  PIV.  With  that  in  mind.  Figure  3-9 
compares  Re  =  20  K,  30  K  and  60  K  dye  injection  in  the  AFRL  water  tunnel.  To  lower  the  operating  Re, 
the  tunnel  is  run  more  slowly,  and  the  rig  and  the  dye  discharge  rate  are  proportionately  slower. 


t/T  =  0 
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t/T  =  0.375 


t/T  =  0.417 


t/T  =  0.75 


Figure  3-9:  AFRL  Water  Tunnel  Dye  Injection  for  Pure-Plunge; 
Re  =  20  K  (left  column),  30  K  (middle  column)  and  60  K  (right). 


At  the  trailing  edge,  at  Re  =  30  K  and  especially  at  Re  =  20  K  there  is  a  discernable  trailing  edge  vortex  at 
t/T  =  0.375  -  0.417,  with  a  region  of  reverse  flow  just  ahead  of  the  trailing  edge.  At  Re  =  60  K,  no  TEV  is 
clearly  visible  in  the  dye  injection,  but  it  was  strongly  apparent  at  t/T  =  0.417  in  the  PIV,  and  at  t/T  =  0.5 
in  the  CFD,  as  mentioned  in  the  previous  section.  The  near-wake  at  the  top  of  the  plunge  stroke  also 
shows  a  discernable  Re-effect,  with  coherent  vortices  seen  for  Re  =  20  K.  In  going  from  t/T  =  0  to  0.125  to 
0.25,  the  suction-side  dye  concentrations  splits,  as  it  were,  into  a  leading-edge  and  trailing-edge  portion, 
the  former  coagulating  into  the  LEV  at  t/T  =  0.25.  Towards  the  bottom  of  the  plunge  stroke,  the  pocket  of 
dye  lacuna  just  aft  of  the  leading  edge  is  smaller  for  higher  Re,  further  suggesting  that  this  region  can  be 
thought  of  as  a  laminar  separation  bubble. 

RANS  computations  and  PIV  measurements  were  performed  at  TU  Darmstadt  at  Re  =  30  K.  In  the  TU 
Darmstadt  wind  tunnel,  rig  max-speed  considerations  precluded  testing  at  Re  =  60  K,  but  Re  =  30  K  was 
readily  attainable.  Vorticity  results  are  presented  in  Figure  3-10  Note  that  in  these  results  the  dimensionless 
vorticity  contour  levels  are  -15  to  +15,  not  -36  to  +36.  Consistent  with  the  dye  injection  at  Re  =  30  K  in 
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Figure  3-9,  there  is  indeed  a  TEV  in  experiments  at  Re  =  30  K,  and  it  is  stronger  in  computation  than  in 
experiment.  But  LEV  formation  does  not  show  much  Re-dependency.  Similar  to  the  UM  RANS,  the  LEV  is 
thinner  and  delayed  in  the  computation  relative  to  the  experiment. 
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Figure  3-10:  TU  Darmstadt  PIV  (left  column)  and  CFD  (right  column)  Vorticity  Contours;  Re  =  30  K; 
Phases  of  Motion  Phi  =  0,  90, 120, 150, 180,  210  and  270;  Contour  Levels  -15  to  +15. 


Turning  to  the  computational  results,  data  are  available  from  the  University  of  Michigan  (Figure  3-11)  and 
the  Middle  East  Technical  University  (Figure  3-12).  For  brevity,  we  focus  on  the  phases  phi  =  90  and  phi  = 
180.  Both  data  sets  are  shown  at  Re  =  60  K,  30  K  and  10  K.  As  Re  decreases,  the  LEV  formation  process  in 
the  UM  computation  begins  to  approach  that  of  the  METU  computation.  Curiously,  and  contradicting  the 
available  experimental  data,  the  TEV  is  stronger  in  UM  computation  at  higher  than  at  lower  Re. 


Figure  3-11:  University  of  Michigan  Computations;  Phases  Phi  =  90  (top  row)  and  180 
(bottom  row);  Re  =  60  K  (left  column),  30  K  (middle  column)  and  10  K  (right  column). 


Figure  3-12:  METU  Computations;  Phases  Phi  =  90  (top  row)  and  180  (bottom  row); 
Re  =  60  K  (left  column),  30  K  (middle  column)  and  10  K  (right  column). 
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Thus  far  we  have  considered  variations  in  the  flowfield  -  resolution  of  LEVs,  “laminar”  vs.  “turbulent” 
separations,  and  so  forth.  Next,  we  turn  to  arguably  the  more  important  question  of  aerodynamic  force 
coefficients  -  lift,  pitch  and  drag.  How  well  to  RANS  methods  predict  force  data,  and  how  much  variation  is 
there  with  Reynolds  number?  We  consider  this  in  the  following  section. 


3.3  AERODYNAMIC  FORCE  COEFFICIENTS 

The  most  basic  quantity  of  interest  is  lift  coefficient,  and  indeed  one  hopes  to  obtain  reasonably  correct  lift 
coefficient  time  history  at  Re  =  60  K  across  the  full  range  of  analytical  and  computational  methods, 
including  the  lower-order  methods.  Figure  3-13  shows  computational  and  measured  lift  coefficient  time 
history  at  Re  =  60  K.  The  encouraging  result  is  that  all  curves  qualitatively  follow  the  same  trend. 
However,  we  note  that  because  of  the  large  peak-to-peak  lift  coefficient  excursion,  relatively  large  errors  - 
say.  Cl  ~  0.4  -  appear  small.  This  is  the  amount  by  which  the  simple  Cl  =  2  Tia  over-predicts  positive  and 
negative  lift  extrema.  Also  interesting  is  that  the  lift  time  history  is  essentially  sinusoidal  and  has  small 
phase  difference  with  respect  to  angle  of  attack  time  history,  despite  the  obvious  presence  of  an  LEV  and 
its  putative  dynamic-stall  effects  on  lift,  which  ought  to  manifest  itself  as  a  large  hysteresis  for  lift  plotted 
vs.  angle  of  attack. 
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Figure  3-13:  Lift  Coefficient  Time  History,  SD7003  Pure-Piunge;  Re  =  60  K; 

Piotted  vs.  Motion  Phase  (ieft)  and  Effective  Angie  of  Attack  (right). 

Experiment  and  AFRL  3D  LES  are  very  close,  as  would  be  expected  from  a  high-fidelity  computation. 
NRC  2D  LES  essentially  splits  the  difference  between  AFRL  LES  and  RANS,  as  is  to  be  expected  from 
the  resolution  of  the  NRC  computation.  But  the  most  remarkable  fact  is  that  Theodorsen’s  formula 
(dashed  black  curve  in  Figure  3-13)  also  follows  fairly  closely  with  the  3D  LES  computations,  slightly 
under-predicting  dynamic  lift  at  the  max  effective  angle  of  attack  (t/T  =  0.25)  and  in  turn  over-predicting 
lift  on  the  first  half  of  the  upstroke  (t/T  =  0.5  to  0.75). 

RANS  computations  predict  slight  lift  coefficient  peaks  near  t/T  ~  0.25,  evidently  due  to  the  LEV.  This  is 
missed  by  the  methods  that  ignore  smooth-surface  flow  separation  -  Theodorsen  and  VPM.  The  closeness 
of  the  VPM  and  Theodorsen  results  implies  that  wake  curvature  is  not  important  for  the  lift  time  history  in 
this  case.  However,  RANS  methods  over-predict  loss  of  lift  at  the  plunge  downstroke  (t/T  =  0.5),  relative 
to  LES  and  experiment,  and  in  general  show  a  slightly  stronger  dynamic  stall  than  does  the  LES;  in  other 
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words,  a  fuller  hysteresis  loop.  Curiously,  the  large  difference  in  velocity/vorticity  contour  plots  between 
the  AFRL  2D  and  3D  LES  and  between  AFRL  and  NRC  LES  does  not  correspond  to  much  difference  in 
the  lift  time  history. 

We  note  that  the  effect  of  LEVs  on  the  lift  time  history  is  modest  -  at  most  10%  -  20%  of  the  maximum 
lift.  This  helps  to  explain  why  lower-order  methods  are  fairly  successful.  But  the  not  unreasonable  result 
from  Cl  =  2  Tia  is  misleading,  because  quasi-steady  thin  airfoil  theory  ignores  unsteady  effects  rather  than 
taking  separation  into  account.  Thus,  whatever  accuracy  to  which  it  may  lay  claim,  this  accuracy  is 
fortuitous.  We  also  note  that  period-to-period  variations  in  the  experiment  may  be  smeared  in  the 
ensemble  average,  in  the  sense  that  dynamic-stall  peaks  and  troughs  in  the  lift  coefficient  time  history  in 
any  one  period  are  attenuated  in  the  average  because  they  vary  randomly  from  period  to  period.  This  is 
disappointing  from  the  viewpoint  of  fundamental  fluid  mechanics,  but  from  the  viewpoint  of  applied 
engineering  one  concludes  that  dynamic  stall  effects  are  uncorrelated  and  therefore  of  limited  importance, 
and  the  lift  time  history  is  quite  sinusoidal,  with  small  phase  lag. 

Variations  in  drag  coefficient  (Figure  3-14)  are  somewhat  larger  than  in  lift.  This  is  to  be  expected,  because 
drag  is  a  smaller  quantity  than  lift  and  therefore  more  sensitive  to  error,  and  because  drag  will  depend  more 
than  lift  on  local  variations  in  pressure  distribution  (presence  or  absence  of  an  LEV,  for  example)  and  on 
prediction  of  boundary  layer  properties.  The  VPM  and  Theodorsen’s  method  are  inviscid,  and  therefore 
disregard  skin  friction  and  LEVs;  but  they  do  capture  Knoller-Betz  type  of  thrust  [27].  Clearly  the  effects  of 
boundary  layers  and  flow  separation  are  evident  in  the  failure  of  Theodorsen’s  formula  (or  more  properly, 
Garrick’s  extension  [27])  to  capture  the  correct  drag  vs.  thrust  balance.  The  UTIAS  method,  which  is 
inviscid  but  has  empirical  corrections  for  viscous  effects,  performs  acceptably  in  phases  of  motion  where  the 
flow  is  nominally  attached,  but  misses  the  effect  of  LEV  formation  and  growth  entirely,  predicting  a  thrust 
near  t/T  -  0.25,  instead  of  a  drag.  AFRL  2D  LES  tends  to  slightly  under-predict  drag  due  to  the  LEV  relative 
to  3D  LES,  and  NRC  2D  LES  under-predicts  further.  The  METU  RANS  computations  track  closely  with  the 
3D  LES,  evidently  because  they  happen  to  resolve  the  LEV  qualitatively  in  a  manner  much  akin  to  the  3D 
LES  vorticity  result.  The  UM  RANS  computation  tracks  closely  with  the  3D  LES  starting  at  t/T  ~  0.4, 
but  predicts  a  slight  thrust  at  the  phase  of  motion  where  the  LEV  forms,  evidently  because  it  does  not 
sufficiently  penalize  the  airfoil  for  loss  of  leading  edge  suction  due  to  leading  edge  separation.  Indeed  all  of 
the  drag  computations,  except  for  UM  RANS,  mutually  track  in  a  narrow  range  except  for  0.25<t/T<0.5, 
which  is  the  period  from  LEV  formation  to  shedding. 


—  - - -  AFRT.  T.K5;  pnie-plniige  2D 


- —  -  .U-'RL  LHS  imic-pluii^c  2D 


Figure  3-14:  Drag  Coefficient  Time  History,  SD7003  Pure-Piunge;  Re  =  60  K; 
Piotted  vs.  Motion  Phase  (ieft)  and  Effective  Angie  of  Attack  (right). 
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Overall  one  can  note  that  whereas  lift  varies  very  close  to  sinusoidally,  despite  all  of  the  various  flow 
separation  effects,  drag  is  not  sinusoidal,  but  rather  peaks  quite  strongly  just  before  the  bottom  of  the 
plunge  downstroke,  where  evidently  flow  separation  is  the  largest.  Drag,  to  belabour  the  obvious,  is  much 
more  sensitive  than  lift  to  the  actual  flowfield  structure,  and  to  calculate  drag  accurately  one  must  resolve 
the  flowfleld  more  accurately  than  is  necessary  for  just  calculating  the  lift. 

We  next  turn  to  aerodynamic  force  trends  with  respect  to  Reynolds  number  (Figure  3-15),  comparing 
AFRL  3D  LES,  UM  RANS  at  Re  =  60  K  -  30  K  -  10  K,  METU  RANS  at  Re  =  60  K  -  30  K  -  10  K, 
and  TUD  RANS  at  Re  =  30  K.  Other  than  the  TUD  lift  relative  overshoot  at  t/T  ~  0.25,  there  is  little 
variation  amongst  the  various  curves,  and  thus  little  lift  dependency  with  Reynolds  number  -  even  less 
than  flowfield  differences  seen  in  for  example  in  Figure  3-9.  A  similar  but  more  tempered  observation 
holds  for  drag.  As  decreases,  the  UM  RANS  drag  prediction  gets  ever  closer  to  the  other  curves  on  the 
second  half  of  the  downstroke,  as  evidently  the  predicted  loss  of  LE  suction  intensifies  with  decreasing 
Reynolds  number.  Low  sensitivity  of  aerodynamic  forces  to  Reynolds  number  is  consistent  with  classical 
dynamic  stall  observations,  for  example  McCroskey  et  al.  [5]. 


.Xf-’RL  LE.^  pme-pliiJitLC 
UM  Re  *  1  OK 

LfM  pme-phiJi^e  Re  =  30K 
UM  pmc-phiJigc  Re  =  5  OK 
piue -plunge  AFRT,  iiie^L^iuieineiil 
METU  piue-pliu^e  Re  =  (jOK 
MET(.’  iJiire-phiJ^ge  Re  ”  3tiK 
METU  ]3iiie-iiliiiige  Re  “  |0K 
TIT>  pme  phuige  Re  3 OK 


.AFRL  LES  piue-pliuige  5D 
LfM  piue-pliiJige  Re  -  l  OR 
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AFRL  LES  pinc-pluiige  3D 
L*M  piue-pliuige  Re  =  JOK 
ITM  pMUC-pliuigc  Re— 30K 
tiM  pue-phuige  Rc  =  t^OK 
l^rETlT  ptiie-phuige  Re  =  fiOK 
^rETU  ]jiiie-|3liiiis'e  Re  =  30K 
METU  pme-iJliuige  Re  =  lOK 
TLTD  pine-pliuis'e  Re  =  30K 


-  AFRL  LE.S  piiie-pliuige  30 

-  UM  jniie-phuige  Re  lOK 

- l^M  piire-phuige  Rc  =  30K 


Figure  3-15:  Trends  in  Lift  Coefficient  (top  row)  and  Drag  Coefficient  (bottom  row) 
for  SD7003  Pure-Piunge  vs.  Reynoids  Number  (60  K,  30  K  and  10  K). 
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The  discrepancy  between  the  UM  and  METU  RANS  codes  deserves  further  mention;  both  are  2D  and  use 
the  same  nominal  turbulence  model.  However,  the  UM  computation  has  two  variations:  “standard”  and 
“modified”,  as  described  in  Chapter  2.  Figure  3-16  compares  the  three  RANS  results  for  lift  and  drag. 
The  “modified”  version  of  the  UM  computation  gives  results  very  close  to  those  to  the  METU  computation. 
While  the  differences  between  all  of  these  curves  are  comparatively  minor  in  the  larger  picture  of  peak-to- 
peak  coefficient  variations,  nevertheless  in  the  details  the  computational  result  is  clearly  sensitive  to  choice 
of  implementation  of  turbulence  model.  It  is  perhaps  disappointing  that  comparatively  significant  variations 
in  aerodynamic  force  prediction  can  arise  from  a  subtlety  of  numerical  implementation. 


Re  —  60K 


Figure  3-16:  Lift  Coefficient  (ieft)  and  Drag  Coefficient  (right)  Comparison 
Between  METU  and  Standard  vs.  Modified  UM  RANS. 

Finally  we  consider  pitching  moment  coefficient.  Cm  is  notoriously  difficult  to  measure  or  compute,  because 
it  depends  not  only  on  the  area  between  the  airfoil  upper-surface  and  lower-surface  pressure  distributions, 
but  on  the  pressure  distributions  directly.  “Correct”  calculation  of  lift  can  sometimes  arise  from  cancellation 
of  errors;  quite  the  opposite  happens  for  pitch,  where  small  errors  can  amplify.  In  experiments,  inertial  tares 
for  Cm  in  a  dynamic  experiment  require  not  only  model  mass  tares,  but  also  moment  of  inertia  tares  -  unless 
testing  in  liquids,  where  inertia  plays  a  much  smaller  role. 

Dynamic  stall  in  classical  applications  such  as  helicopter  blades  and  manoeuvring  aircraft  is  even  more 
notorious  for  its  effect  on  pitching  moment,  than  on  lift.  Formation  of  the  LEV  in  dynamic  stall  causes  a 
small  nose-down  pitching  moment  relative  to  the  quarter-chord,  because  the  LEV  is  close  to  the  quarter- 
chord  and  its  local  suction  has  a  small  moment  arm.  However,  as  the  LEV  convects  downstream  and  reaches 
the  trailing  edge,  the  pitching  moment  coefficient  becomes  strongly  nose-up.  Then,  when  the  vortex  is  aft  of 
the  airfoil.  Cm  snaps  back  up,  producing  a  violent  hysteresis  when  plotted  against  a  [6].  It  is  also  possible 
that  the  strong  trailing  edge  vortex  found  in  the  computations,  if  present  in  the  flow,  would  add  to  a  moment- 
stall  as  it  forms  and  convects  downstream  from  the  trailing  edge. 

The  available  Cm  data  are  shown  in  Figure  3-17.  Clearly  there  is  a  classical  dynamic  moment  stall  - 
in  contradistinction  to  the  situation  for  lift.  This  is  missed  by  the  inviscid  methods,  because  these  methods 
like  a  LEV  model.  Because  the  METU  RANS  have  qualitatively  similar  LEV  resolution  to  AFRL  3D  EES, 
the  two  Cm  results  are  similar.  2D  AFRL  EES  does  not  predict  a  coherent  LEV,  whence  there  is  much  less 
moment  stall.  Cm  prediction,  completely  unlike  Cl  and  to  a  large  extent  unlike  Cd,  requires  accurate 
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rendition  of  the  flowfield.  But  that  accuracy,  if  achieved  by  a  low-order  method,  can  occasionally  be 
serendipitous.  We  mention  the  trite  but  true  platitude  that  test  cases  across  a  broad  range  of  motion 
parameters  are  necessary,  before  one  can  speak  definitively  about  the  accuracy  of  low-order  methods. 


Figure  3-17:  Pitching  Moment  Coefficient  for  SD7003  Airfoii  in  Pure-Piunge. 


To  summarize,  an  effective  peak  angle  of  attack  of  22°  at  Re  =  60  K  is  large  enough  to  cause  moment-stall  in 
the  sense  of  classical  dynamic  stall,  but  not  lift  stall.  Boundary  layer  physics  does  not  appear  to  be  a  main 
driver  of  integrated  aerodynamic  effects.  This  is  largely  positive  news  for  engineering-level  predictive 
methods.  Clearly  a  larger  parameter  study  is  warranted,  exploring  for  example  the  role  of  airfoil  LE  radius 
and  its  effect  on  boundary  layer  physics,  and  looking  at  more  aggressive  motions  with  larger  effective  angle 
of  attack.  Presently,  however,  we  turn  to  the  arguably  more  complicated  problem  of  combined  pitch-plunge, 
where  the  effective  angle  of  attack  is  smaller  due  to  partial  cancellation  of  plunge-induced  angle  of  attack 
from  pitch.  As  we  will  see,  the  resulting  flow  physics  are  rather  less  spectacular,  but  perhaps  more  complex 
from  the  viewpoint  of  engineering  prediction. 
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Chapter  4  -  COMBINED  PITCH-PLUNGE 
OF  THE  SD7003  AIRFOIL 


In  contrast  with  the  pure-plunge  motion,  the  combined  pitch-plunge  motion  is  at  effective  angle  of  attack 
not  far  above  static  stall,  with  the  motivation  of  serving  as  a  case  for  useful  propulsive  thrust. 


4.1  RESUME  OF  EXPERIMENTAL  AND  COMPUTATIONAL  RESULTS  FOR 
THE  FLOWFIELD:  VELOCITY  AND  VORTICITY 

As  for  the  pure-plunge  case,  we  first  consider  velocity  and  vorticity  contours,  and  then  proceed  to  comparing 
results  for  aerodynamic  loads.  The  canonical  Reynolds  number  remains  at  60000. 

4.1.1  Computations 

RANS  computations  were  conducted  by  METU  (Figure  4-1)  and  UM  (Figure  4-3).  NRC  performed  2D  LES 
(Figure  4-2). 
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Figure  4-1:  METU  Computations;  Phases  Phi  =  0,  90, 120, 150, 180,  210  and  270. 
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Figure  4-2:  lAR  NRC  Computations,  LES;  Phases  Phi  =  0,  90, 180  and  270. 
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Figure  4-3:  University  of  Michigan  Computations;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 


4.1.2  Experiments 

Experiments  were  performed  by  UM  (Figure  4-4),  AFRL  (three  iterations  spaced  over  two  years:  Figure  4-5, 
Figure  4-6  and  Figure  4-7,  respectively),  and  TUBS  (Figure  4-8). 
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Figure  4-4:  University  of  Michigan  PiV;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 


Figure  4-5:  AFRL  PiV,  First  Data  Series;  Phases  Phi  =  0,  90, 180  and  270. 
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Figure  4-6:  AFRL  PIV,  Second  Data  Series;  Phases  Phi  =  0,  45,  90, 135, 180,  225,  270  and  315. 
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Figure  4-7:  AFRL  PIV,  Third  Data  Series;  Phases  Phi  =  0,  90, 120, 150, 180,  210  and  270. 


Figure  4-8:  TUBS  PiV;  Phases  Phi  =  0,  90,  180  and  270. 
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Reviewing  the  streamwise  velocity  contour  plots  for  the  pitch-plunge  case,  results  group  into  three  types: 
“small”  separation,  where  the  region  of  stagnant  or  reversed  flow  at  the  bottom  of  the  plunge  downstroke  is 
thin  and  does  not  begin  until  towards  the  airfoil  mid-chord;  “large”  separation,  which  is  still  closed  in  the 
phase-averaged  sense  but  begins  smoothly  shortly  aft  of  the  leading  edge,  and  LEV-dominated  large 
separation.  The  U  of  Michigan  experiment  and  computation,  and  one  of  the  AFRL  experiments  fall  under 
“small”  separation.  The  NRC  computation,  TUBS  experiment  and  another  AFRL  experiment  fall  under  the 
“large”  separation.  The  METU  computation  shows  an  LEV. 

The  disparity  between  the  three  kinds  of  separation  is  loosely  associable  with  laminar  to  turbulent  transition. 
The  METU  computations  appear  to  me  “more  laminar”,  while  the  U  of  Michigan  computations  are  “more 
turbulent”.  The  AFRL  water  tunnel  is  evidently  on  the  edge  between  “laminar”  and  “turbulent”  results, 
sometimes  favouring  the  one  and  on  occasion  the  other. 

It  is  perhaps  worthwhile  to  briefly  note  the  factors  behind  turbulent-like  or  laminar-like  results.  In  the 
computations,  causes  include  of  course  the  choice  of  turbulence  model,  the  inputs  to  that  model,  the  role 
of  the  mesh  and  treatment  of  far-field  boundary  conditions,  the  relative  balance  between  turbulent  kinetic 
energy  production  and  dissipation,  and  so  forth.  For  experiments,  the  main  cause  is  disturbances,  which 
can  be  flowfield  related  or  model-related.  Flowfield  disturbances  include  spatial  variation  of  flow  speed 
and  possible  secondary  flows  in  the  test  section,  temporal  variation  (“turbulence  intensity”),  blockage 
effects  and  streamline  curvature  and  change  in  the  effective  angle  of  attack,  vortical  disturbances  that 
survive  the  screens  in  the  settling  chamber  and  enter  the  test  section,  and  on  and  on.  Model-related 
problems  are  especially  acute  for  dynamic  testing.  The  include  elastic  deflections  and  vibrations  of  the 
model  or  the  rig  linkages,  departure  of  actual  model  motion  from  the  prescribed  motion,  vibrations,  model 
surface  imperfections  and  so  forth.  The  list  is  vast  and  ambiguous,  as  indeed  the  recitation  of  problems 
basically  speaks  to  aerodynamic  computation  and  experiment  in  general,  not  just  at  low  Reynolds  number 
and  not  just  for  dynamic  testing.  The  point  here  is  not  to  attempt  to  identify  flowfield  features  with 
specific  factors  such  as  drawbacks  of  some  turbulence  model  or  blockage  in  a  wind  tunnel,  but  rather  to 
note  that  the  present  case,  where  the  effective  angle  of  attack  just  slightly  exceeds  the  static  stall  angle  of 
the  SD7003  airfoil,  is  very  sensitive  to  “transition”  in  general.  If  one  wishes  to  optimize  for  flapping-wing 
cruising  propulsion,  say,  and  it  is  important  to  capture  the  aerodynamic  coefficient  time  history 
quantitatively,  then  RANS-type  computations  and  ground- facility  tests  would  alike  require  great  care. 
This  is  in  marked  contrast  with  the  deep-stall  case  in  the  previous  chapter. 

However,  the  meaning  of  turbulence  can  be  subtle;  smaller  free-stream  turbulence  leads  to  larger  separation, 
and  inside  the  larger  separated  flowfield  structure  the  turbulent  kinetic  energy  can  be  larger  than  in  the 
smaller  one.  As  a  rough  comparison,  contours  of  2D  turbulent  kinetic  energy  from  two  AFRL  PIV  data  sets, 
both  at  the  bottom  of  the  plunge  stroke,  are  given  in  Figure  4-9. 
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Figure  4-9:  PIV-Derived  Planar  Turbulent  Kinetic  Energy  Contours,  AFRL  Data  Sets;  Phase  Phi  =  180 
(Bottom  of  Plunge  Downstroke);  “Small”  Separation  (left)  and  “Large”  Separation  (right). 
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Another  potential  area  of  concern  is  spanwise  variations.  The  spanwise  extent  of  flow  separation  for  both  the 
pure-plunge  and  pitch-plunge  SD7003  cases  was  reported  in  01  et  al.  [55],  and  indeed  this  has  implications 
for  the  role  of  blockage  (for  experiments)  and  spanwise  extent  of  the  domain  (for  computations).  But  a 
related  question,  not  explored  in  01  et  al.  [55],  is  the  extent  to  which  the  flow  separation  varies  in  going 
across  the  model  span  towards  the  tunnel  wall;  that  is,  how  prominent  are  the  sidewall  effects?  If  the  PIV 
light  sheet  is  placed  too  close  to  the  tunnel  sidewalls,  then  possibly  the  measured  flowfleld  will  depart  from 
that  in  the  “core”  of  the  test  section  flow.  Some  evidence  is  available  from  the  cross-flow-plane  PIV  data  of 
TU  Braunschweig,  where  contours  of  streamwise  velocity  and  normalized  Reynolds  stress  are  plotted  in 
Figure  4-10,  for  the  chordwise  location  of  x/c  =  0.40,  at  the  bottom  of  the  plunge  downstroke.  Data  are  taken 
from  nearly  at  the  tunnel  test  section  starboard  wall,  towards  -0.45  of  the  airfoil  span.  It  is  evident  that  the 
region  around  0.2  or  0.25  spans  away  from  the  wall  is  a  dividing  point;  inboard  of  this  region,  the  flow  is 
spanwise-uniform,  while  outboard  of  this  region  and  approach  the  tunnel  wall,  there  is  a  spanwise  variation. 
The  implication  is  important  for  example  for  the  AFRL  experiments;  since  there  the  PIV  light  sheet  was 
approximately  at  the  %  span  location,  slight  variations  further  inboard  or  outboard  could  materially  affect  the 
resulting  data.  This,  presumably,  is  even  more  true  for  the  more  separated  flow  in  the  pure-plunge  case, 
but  unfortunately  data  analogous  to  Figure  4-10  is  not  available  for  that  case. 
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Figure  4-10:  TU  Braunschweig  PIV;  Contours  of  Streamwise  Velocity  (left)  and  Normalized 
Reynolds  Stress,  ~u’v’  (right),  at  Motion  Phase  of  180  Degrees  (Bottom  of 
Plunge  Downstroke);  Chordwise  Location  x/c  =  0.40. 


4.2  REYNOLDS  NUMBER  EFFECTS 

As  for  the  pure-plunge  case  presented  in  the  previous  chapter,  we  consider  Reynolds  numbers  below  60  K. 
Such  an  assessment  is  even  more  important  for  the  lower  total  effective  angle  of  attack  case  considered 
presently,  because  presumably  natural  transition  effects  should  be  less  overwhelmed  by  high  pressure 
gradients  at  the  leading  edge,  caused  by  the  airfoil  motion.  At  Re  =  30  K  we  again  consider  vorticity 
contour  plots  for  the  Darmstadt  PIV  vs.  RANS  comparison  (Figure  4-1 1),  with  contour  levels  -15  to  +15. 
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Figure  4-11:  TU  Darmstadt  PIV  (left  column)  and  CFD  (right  column)  Vorticity  Contours;  Phases 
of  Motion  Phi  =  0,  90, 120, 150, 180,  210  and  270;  Re  =  30  K;  Contour  Levels  -15  to  +15. 
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Next  we  compare  the  series  of  RANS  computations  from  the  University  of  Michigan  (Figure  4-12)  and 
Middle  East  Technical  University  (Figure  4-13),  for  phases  phi  =  90  and  180,  at  Re  =  60  K,  30  K  and 
10  K.  In  all  cases  there  is  a  general  trend  towards  larger  flow  separation  at  lower  Re.  This  is  confirmed  by 
the  AFRL  water  tunnel  dye  injection  results,  which  were  timed  to  find  the  “small”  separation  at  Re  =  60  K 
(Figure  4-14). 


Figure  4-12:  University  of  Michigan  Computations;  Phases  Phi  =  90  (top  row)  and  180 
(bottom  row);  Re  =  60  K  (ieft  coiumn),  30  K  (middie  coiumn)  and  10  K  (right  coiumn). 


Figure  4-13:  Middie  East  Technicai  University  Computations;  Phases  Phi  =  90  (top  row)  and 
180  (bottom  row);  Re  =  60  K  (ieft  coiumn),  30  K  (middie  coiumn)  and  10  K  (right  coiumn). 


t/T  =  0 
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t/T  =  0.75 


Figure  4-14:  AFRL  Water  Tunnel  Dye  Injection;  Phases  Phi  =  0,  90, 120, 180  and  270; 

Re  =  10  K  (left  column),  30  K  (middle  column)  and  60  K  (right  column). 

We  next  consider  how  the  various  differences  in  the  flowfield  -  due  to  computational  methods  or 
experimental  setups,  and  Reynolds  number  -  translate  into  differences  in  aerodynamic  force  coefficients. 


4.3  AERODYNAMIC  FORCE  COEFFICIENTS 

As  with  the  pure-plunging  case,  we  are  interested  in  comparing  the  flowfield  evolution  history  with  the 
aerodynamic  coefficient  time  history,  again  principally  the  lift.  Lift  coefficient,  plotted  vs.  motion  phase  and 
vs.  effective  angle  of  attack,  is  given  in  Figure  4-15,  while  drag  is  shown  in  Figure  4-16.  Once  again, 
the  quasi-steady  approximation  and  Theodorsen’s  formula  are  the  simplest  models,  followed  by  the  vortex- 
particle  method,  then  by  the  two  RANS  computations,  from  the  University  of  Michigan  and  METU, 
and  finally  by  the  LES  computation  (NRC). 
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Figure  4-15:  Lift  Coefficient  Time  History,  SD7003  Pitch-Plunge;  Re  =  60  K; 
Plotted  vs.  Motion  Phase  (left)  and  Effective  Angle  of  Attack  (right). 
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Figure  4-16:  Drag  Coefficient  Time  History,  SD7003  Pitch-Piunge;  Re  =  60  K; 

Piotted  vs.  Motion  Phase  (ieft)  and  Effective  Angie  of  Attack  (right). 

With  this  large  of  a  spread  in  data,  it  is  difficult  to  make  a  statement  about  method  complexity  vs.  method 
accuracy,  to  the  level  possible  with  the  pure -plunge  case.  Nevertheless,  one  can  make  some  qualitative 
observations.  Comparing  the  two  RANS  computations,  we  find  that  indeed  the  more  “laminar  like” 
(METU)  predicts  stronger  vortical  effects  on  ClO),  than  does  the  more  “turbulent  like”  (UM).  The  UM 
RANS  computation  and  the  NRC  LES  are  quite  close,  and  the  AFRL  experiment  is  not  far  from  the  two. 
Theodorsen  over-predicts  lift  in  t/T<0.25<0.75,  and  interestingly  enough,  the  Theodorsen  prediction 
appears  worse  for  pitch-plunge  than  for  pure-plunge! 

The  UTIAS  lumped-vortex  and  AFRL  VPM  are  conceptually  close,  but  the  UTIAS  result  tracks  closer 
with  what  might  be  termed  the  average  lift  time  history,  evidently  because  it  has  a  static  and  dynamic  stall 
model,  whereas  the  AFRL  vortex  particle  method  is  a  classical  attached-flow  approach.  This  probably  also 
explains  why  the  UTIAS  method  appears  to  perform  well  in  drag  prediction  (Figure  4-16);  it  tracks  very 
closely  with  the  UM  RANS,  which  of  course  is  conceptually  an  entirely  different  approach.  The  NRC  2D 
LES  predicts  a  lower  thrust  at  t/T  ~  0.25,  evidently  due  to  the  larger  predicted  separation.  Contrary  to  the 
pure-plunge  case  discussed  in  the  previous  chapter,  Theodorsen’ s  formula  is  fairly  successful  in  at  least 
qualitative  prediction  of  the  drag  coefficient  history,  except  late  in  the  motion  upstroke,  where  evidently 
turbulent  skin  friction  drag  is  main  cause  of  discrepancy. 

We  next  consider  dependency  of  lift  and  drag  with  Reynolds  number  (Figure  4-17),  again  comparing  three 
RANS  computations  (UM,  METU  and  TUD).  Unlike  the  pure-plunge  case,  reduction  in  Re  from  60  K  to 
30  K  to  10  K  does  show  significant  changes  in  lift  and  drag  and  significant  scatter  between  the  computations. 
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Figure  4-17:  Trends  in  Lift  Coefficient  (top  row)  and  Drag  Coefficient  (bottom  row) 
for  SD7003  Pitch-Piunge  vs.  Reynoids  Number  (60  K,  30  K  and  10  K). 


As  with  the  pure -plunge  case,  we  compare  the  “standard”  and  “modified”  UM  RANS  with  the  METU 
RANS,  to  assess  the  significance  of  implementation  of  nominally  the  same  turbulence  model.  Again, 
the  “modified”  version  of  the  UM  computation  agrees  very  closely  with  the  METU  computation. 
But  unlike  the  pure-plunge  case,  the  respective  differences  are  qualitatively  significant,  and  are  not  merely 
wiggles  here  and  there.  This,  in  the  shallow-stall  case  even  more  caution  is  required  in  turbulence  model 
implementation  than  in  the  deep-stall  case.  As  we  shall  see  in  the  next  chapter,  the  differences  between  the 
various  results  are  significantly  attenuated  when  there  is  a  forcing  of  separation  (and  presumably 
transition)  due  to  sectional  geometry,  as  we  turn  to  the  case  of  a  fiat  plate  of  round  leading  edges. 
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Figure  4-18:  Lift  Coefficient  (ieft)  and  Drag  Coefficient  (right)  Comparison 
Between  METU  and  Standard  vs.  Modified  UM  RANS. 

The  pitch-plunge  case  was  selected  not  only  to  explore  effects  of  Reynolds  number,  transition,  turbulence 
models,  facility  flow  quality  and  so  forth  -  but  to  consider  a  case  suitable  for  propulsive  flapping-wing  flight, 
in  the  sense  of  flight  of  large  birds,  as  discussed  in  Chapter  1.  Does  the  subject  motion  indeed  produce  net 
thrust,  where  Garrick-type  propulsion  from  leading  edge  suction  [27]  overcomes  friction  and  pressure 
(laminar  separation  bubble,  separated  wake)  drag?  The  flndings  are  summarized  in  Table  4-1.  Unfortunately 
there  are  no  experimental  results  for  drag,  but  all  of  computations  that  solve  the  Navier-Stokes  equations 
predict  a  net  drag.  The  UM  RANS  at  Re  =  60  K  comes  close  to  net  zero  drag,  and  possibly  at  slightly  higher 
Reynolds  number  there  would  be  a  net  thrust.  It  appears,  therefore,  that  even  at  Re  =  60  K  the  motion  studied 
here  is  not  net-propulsive  because  viscous  effects  overwhelm  leading  edge  suction. 

Table  4-1:  Mean  Lift  and  Drag  Coefficients  for  SD7003  Pitch-Piunge 


Data  Set 

Mean  Cl 

Mean  C] 

UMRANSRe=  lOK 

0.694 

0.032 

UM  RANS  Re  =  30  K 

0.839 

0.011 

UM  RANS  Re  =  60  K 

0.894 

0.004 

METU  RANS  Re  =  60  K 

0.856 

0.022 

METU  RANS  Re  =  30  K 

0.801 

0.038 

METU  RANS  Re  =  10  K 

0.682 

0.061 

NRC  EES 

0.849 

0.023 

UTIAS 

0.882 

-0.008 

TUD  RANS  Re  =  30  K 

0.821 

0.028 

Theodorsen 

0.985 

-0.046 
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Chapter  5  -  PLUNGE  AND  PITCH-PLUNGE 
OF  A  FLAT  PLATE 


Here  our  objective  is  to  assess  the  role  of  sectional  geometry,  by  comparing  a  flat  plate  with  the  SD7003 
airfoil.  The  canonical  Reynolds  number  remains  at  60000,  and  the  kinematics  of  motion  are  the  same  as  in 
the  prior  two  chapters. 


5.1  RESUME  OF  EXPERIMENTAL  AND  COMPUTATIONAL  RESULTS  FOR 
THE  FLOWFIELD:  VELOCITY  AND  VORTICITY 

5.1.1  Computations,  Pitching-Plunging  Plate 

RANS  computations  were  run  by  the  University  of  Michigan  (Figure  5-1)  and  Middle  East  Technical 
University  (Figure  5-2).  Clearly  there  is  less  difference  in  the  velocity  and  vorticity  contours  for  these  two 
computations  for  the  flat-plate,  than  there  were  for  the  corresponding  SD7003  case. 
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Figure  5-1:  University  of  Michigan  Computations;  Phases  Phi  =  0,  90,  120, 150, 180  and  270. 
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Figure  5-2:  Flat-Plate  Pitch-Plunge,  METU  Computations,  Fluent, 
2D;  Phases  Phi  =  0,  90,  120, 150,  180,  210  and  270. 


5.1.2  Experiments,  Pitching-Plunging  Plate 

Experiments  were  limited  to  water  tunnels:  AFRL  (Figure  5-3)  and  University  of  Michigan  (Figure  5-4). 
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Figure  5-3:  AFRL  Water  Tunnel;  Phases  Phi  =  0,  90, 180  and  270. 


Figure  5-4:  University  of  Michigan  Water  Tunnel;  Phases  Phi  =  0,  90, 180  and  270. 
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5.2  REYNOLDS  NUMBER  EFFECTS 

The  flat-plate  is  much  less  sensitive  to  Reynolds  number  than  is  the  SD7003  airfoil,  because  the  plate’s 
leading  edge  essentially  forces  transition,  and  because  the  plate  has  a  very  narrow  range  of  angle  of  attack 
over  which  flow  remains  fully  attached  in  the  static  case.  This  assertion  is  essentially  confirmed  in  the  dye 
injection  survey  at  Re  =  20  K  and  60  K  in  Figure  5-5. 


Figure  5-5:  AFRL  Water  Tunnel  Dye  Injection  for  Pitching-Plunging  Flat  Plate;  Phases  Phi  =  0, 
45,  90,  120,  150,  180,  210  and  270;  Re  =  20  K  (left  column)  and  Re  =  60  K  (right  column). 
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5.2.1  Computations,  Plate  in  Pure-Plunge 

Computations  for  the  wall-to-wall  plate  in  pure-plunge  were  pursued  by  METU  (Figure  5-6)  and  UM 
(Figure  5-7).  The  similarity  between  the  two  is  quite  close,  and  is  much  closer  than  for  the  SD7003  airfoil 
executing  the  same  motion  at  the  same  Re. 
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Figure  5-6:  Flat-Plate  Plunge,  METU  Computations,  Fluent, 
2D;  Phases  Phi  =  0,  90,  120,  150,  180,  210  and  270. 
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Figure  5-7:  University  of  Michigan  Computations,  Waii-to-Waii 
Piunging  Fiat  Piate;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 


5.2.2  Experiments,  Plate  in  Pure-Plunge 

Unfortunately  only  one  series  of  experiments  for  the  wall-to-wall  plate  in  pure-plunge  was  available  - 
from  the  University  of  Michigan  water  tunnel  (Figure  5-8).  As  compared  to  experiments  in  the  same 
facility  for  the  SD7003  airfoil  in  pure-plunge,  the  differences  are  not  great;  the  flat-plate  LEV  is  somewhat 
larger  and  more  diffuse  than  for  the  airfoil,  and  the  low-speed  region  on  the  suction  side  towards  the 
bottom  of  the  downstroke  persists  less. 


Figure  5-8:  University  of  Michigan  PiV,  Waii-to-Waii  Piunging 
Fiat  Piate;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 
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5.3  AERODYNAMIC  FORCE  COEFFICIENTS  FOR  THE  FLAT  PLATE 

5.3.1  Pure-Plunge 

Following  experience  with  the  SD7003  airfoil  in  pure-plunge,  we  expect  at  most  moderate  variation  in 
force  coefficient  time  history  amongst  the  various  methods,  at  least  for  lift.  Indeed,  in  Figure  5-9  the  two 
RANS  computations  are  very  close.  No  LES  for  the  flat-plate  in  pure-plunge  was  available,  but  simply 
reprinting  the  AFRL  3D  LES  result  for  the  plunging  SD7003  shows  close  tracking  with  the  two  plate 
RANS  results.  It  is  possible,  in  fact,  the  RANS  is  over-predicting  the  dynamic  stall  peak  (t/T  ~  0.25). 
That  said,  the  quasi-steady  result  is  quite  accurate  on  the  downstroke,  while  Theodorsen  is  more  accurate 
on  the  upstroke,  respectively. 
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Figure  5-9:  Lift  Coefficient  History  for  the  Flat  Plate  in  Pure-Plunge,  Plotted  vs.  Time  (left) 
and  Angle  of  Attack  (right).  Theodorsen  and  quasi-steady  calculations  have  camber 
removed  and  the  AFRL  LES  is  a  reprint  of  the  respective  SD7003  result. 


Trends  for  drag  are  shown  Figure  5-10.  Again  the  two  RANS  computations  are  very  close,  and  both  are 
close  to  the  SD7003  LES  from  t/T  ~  0.3  onward. 
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Figure  5-10:  Drag  Coefficient  History  for  the  Fiat  Piate  in  Pure-Piunge,  Piotted  vs.  Time  (ieft) 
and  Angie  of  Attack  (right)  -  The  AFRL  LES  is  a  Reprint  of  the  Respective  SD7003  Resuit. 


5.3.2  Pitch-Plunge 

Finally  we  consider  aerodynamic  force  coefficients  for  the  wall-to-wall  flat  plate  in  pitch-plunge;  that  is, 
X  =  0.6:  lift  (Figure  5-11)  and  drag  (Figure  5-12).  For  this  case,  one  LES  computation  is  available, 
from  NRC,  in  2D.  Interestingly  enough,  the  LES  tracks  very  closely  with  the  vortex  particle  method  solution 
and  Theodorsen.  The  quasi-steady  solution,  as  usual,  over-estimates  lift  at  both  crest  (t/T  =  0.25)  and  trough 
(t/T  =  0.75).  The  two  RANS  computations,  METU  and  UM,  are  very  close,  and  predict  much  stronger  lift 
peak  due  to  the  LEV  and  loss  of  lift  as  the  LEV  sweeps  downstream  away  from  the  plate  suction  side. 


AFRL  VPM  pitcJi  pliuifie  no  CLunbei 


Figure  5-11:  Lift  Coefficient  History  for  the  Fiat  Piate  in  Pitch-Piunge,  Piotted  vs.  Time  (ieft)  and 
Angie  of  Attack  (right)  -  Theodorsen  and  Quasi-Steady  Caicuiations  have  Camber  Removed. 
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Figure  5-12:  Drag  Coefficient  History  for  the  Fiat  Piate  in  Pitch-Piunge, 
Piotted  vs.  Time  (ieft)  and  Angie  of  Attack  (right). 


In  drag,  we  are  as  usual  left  with  the  viscous-only  methods,  and  unfortunately  experimental  data  is  not 
available.  The  two  RANS  results  are  again  very  close,  and  only  deviate  from  the  LES  significantly  at  the 
start  of  the  upstroke.  As  compared  with  the  pitching-plunging  airfoil,  the  drag  trough  at  t/T  -  0.25  is  much 
less  acute  for  the  flat-plate,  consistent  with  the  separation  on  the  suction  side  being  more  severe. 

Finally,  it  is  worth  mentioning  a  word  on  the  RANS  solutions  for  the  flat  plate,  especially  since  RANS  is  the 
most  common  computational  scheme  for  these  sorts  of  problems.  Whatever  the  differences  of  turbulence 
modelling,  production  and  dissipation  between  the  METU  and  UM  RANS  approaches,  these  differences  are 
essentially  moot  for  the  flat  plate  with  round  leading  edge,  evidently  because  this  leading  edge  is  strongly 
forcing  towards  transition.  That  is,  it  may  be  premature  to  claim  with  treatment  of  turbulence  is  “better”  in 
which  code,  and  in  which  circumstance  -  but  one  can  observe  that  whatever  the  various  differences  may  be, 
the  upshot  is  essentially  no  difference  for  strongly  forced  flows.  It  is  likely,  further,  that  if  we  were  to  run  a 
parameter  study  on  k  and  kh,  and  looked  at  more  violent  motions,  then  the  differences  between  the  various 
approaches  would  be  even  smaller.  This  idea  is  pursued,  if  only  episodically,  in  a  subsequent  chapter. 
But  first  we  turn  to  finite-span  effects,  retaining  the  same  sectional  shape  (thin  flat  plate  with  round  edges), 
but  now  at  AR  =  2. 
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Chapter  6  -  FINITE  ASPECT  RATIO  EFFECTS:  PLUNGE  AND 
PITCH-PLUNGE  OF  AN  ASPECT  RATIO  2  FLAT  PLATE 

6.1  MOTIVATIONS  AND  PROBLEM  DEFINITION 

Perhaps  the  most  natural  generalization  from  wall-to-wall  models  (in  experimental  facilities)  and  spanwise- 
periodic  boundary  conditions  (in  computations)  is  the  finite-wing  of  high  aspect  ratio.  By  systematically 
reducing  aspect  ratio  from  some  large  but  finite  value,  towards  ever  smaller  values,  it  is  possible  to  construct 
a  consistent  passage  into  3D  from  2D.  This  is  sensible,  for  example,  in  terms  of  understanding  the  flapping- 
flight  of  large  birds  (see  Chapter  1),  which  tend  to  be  of  high  aspect  ratio.  However,  in  an  experimental 
facility  this  introduces  the  problem  of  blockage  or  very  low  Reynolds  number.  To  avoid  blockage, 
the  wingspan  needs  to  be  some  fraction  of  the  test  section  width,  and  at  high  aspect  ratio  this  severely  limits 
the  chord.  In  turn  the  Reynolds  number,  based  on  chord,  becomes  very  low.  But  very  low  Reynolds  number 
problems,  besides  being  inconsistent  with  the  focus  of  earlier  chapters,  tend  to  be  of  comparatively  low 
aspect  ratio,  as  is  the  case  for  most  insects.  In  computation,  full  resolution  of  high  aspect  ratio  means  a  grid 
much  longer  in  the  spanwise  than  in  the  streamwise  direction,  which  raises  problems  of  computational  size, 
if  we  wish  to  maintain  high  resolution  in  the  streamwise  direction. 

The  alternative  is  to  consider  low  aspect  ratios.  This  is  more  amenable  to  investigation,  but  more 
importantly,  it  is  especially  interesting  because  of  the  strong  interaction  expected  at  low  aspect  ratio  between 
leading  edge  vortices  and  wingtip  vortices,  where  for  example  the  latter  might  stabilize  the  former  through 
spanwise  pressure  gradients.  Low  aspect  ratio  wings  undergoing  various  longitudinal-plane  manoeuvres  are 
expected  to  evince  a  qualitatively  different  flow  structure  and  loads  time  history  from  those  of  the  same 
sectional  geometry,  but  in  2D.  A  full  treatment  of  the  problem  requires  detailed  3D  velocimetry,  which 
unfortunately  is  beyond  the  scope  of  this  study.  In  keeping  with  the  above-reported  results,  flowfield  data  are 
presented  in  streamwise-parallel  planes,  typically  at  the  %-span  location  of  the  model.  The  featured 
configuration  is  the  rectangular-planform  fiat  plate  of  aspect  ratio  2,  with  round  edges  on  its  entire  periphery. 
Both  the  pure-plunge  and  the  pitch-plunge  case  are  considered.  The  canonical  Reynolds  number  is  now 
40000  rather  than  the  60000  for  the  wall-to-wall  cases,  as  blockage  problems  are  reduced  with  a  smaller 
model,  and  thus  a  smaller  chord. 

The  sectional  geometry  of  a  thin  flat  plate  with  round  edges  was  selected  to  deemphasize  the  role  of 
boundary  layer  transition  and  to  at  least  partially  solve  the  problem  of  how  to  treat  the  airfoil  section  at  the 
wingtips.  Such  shapes  are  also  easy  to  manufacture,  and  in  general  easy  to  grid.  Looking  towards  future 
work,  they  are  amenable  to  generalization  to  flexible  structures,  such  as  membranes. 

We  begin  with  the  pure-plunge  kinematics,  again  considering  velocity  and  vorticity  contour  plots,  and  load 
time  history;  then  we  continue  with  pitch-plunge. 


6.2  THE  PURE-PLUNGE  CASE 

6.2.1  Computations,  AR  =  2  Plate  in  Pure-Plunge 
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Figure  6-1:  U  of  Michigan  Computations,  AR  =  2  Piate  in 
Pure-Piunge;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 


6.2.2  Experiments,  AR  =  2  Plate  in  Pure-Plunge 

Experiments  were  run  in  the  UM  and  AFRL  water  tunnels  (Figure  6-2  and  Figure  6-3,  respectively),  and  in 
the  U  of  Florida  REEF  wing  tunnel  (Figure  6-4). 
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Figure  6-2:  PIV,  U  of  Michigan,  AR  =  2  Piate  in  Pure-Piunge;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 
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Figure  6-3:  PIV,  AFRL,  AR  =  2  Plate  in  Pure-Plunge;  Phases  Phi  =  0,  90, 120, 150, 180,  210  and  270. 


Figure  6-4:  PIV,  U  of  Florida  REEF,  AR  =  2  Plate  in  Pure-Plunge;  Phases  Phi  =  0,  90,  180  and  270. 
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ORGAISriZATION 


6.3  THE  PITCH-PLUNGE  CASE 

6.3.1  Computations,  Pitching-Plunging  Plate,  AR  =  2 


Figure  6-5:  U  of  Michigan  Computations,  Pitching-Piunging 
AR  =  2  Piate;  Phases  Phi  =  0,  90, 120, 150, 180  and  270. 
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6.3.2  Experiments,  Pitching-Plunging  Plate,  AR  =  2 

Flowfield  velocity  and  vorticity  contours  (PIV)  were  obtained  from  the  University  of  Florida  (Figure  6-6), 
AFRL  (Figure  6-7),  and  the  University  of  Michigan  (Figure  6-8). 


Figure  6-6:  U  of  F  REEF  PIV  Results  for  AR  =  2  Plate  Pitch-Plunge. 
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Figure  6-7:  AFRL  PIV  Results  for  AR  =  2  Plate  Pitch-Plunge. 
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Figure  6-8:  U  of  Michigan  PiV  Resuits  for  AR  =  2  Piate  Pitch-Piunge. 


6.4  AERODYNAMIC  FORCE  TIME  HISTORIES  FOR  THE  AR  =  2  FLAT 
PLATE 

For  AR  =  2,  both  slender  body  theory  and  lifting  line  theory  give  a  quasi-steady  lift  curve  slope  of  n.  For  the 
steady  case  at  angles  of  attack  below  stall,  Kaplan  et  al.  [53]  measured  a  lift  curve  slope  for  a  square-edged 
AR  =  2  plate  very  close  to  tt,  with  almost  negligible  non-linear  deviations  due  to  vortical  effects  such  as  loss 
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of  leading-edge  suction.  The  principal  question  now  is  whether  the  AR  =  2  plate’s  behaviour  in  the  unsteady 
problem  is  similarly  simple. 

6.4.1  Forces  for  the  AR  =  2  Plate  in  Pitch-Plunge 

Lift  coefficient  results  for  the  pitching-plunging  case  are  shown  in  Figure  6-9.  Force  measurements  from 
TUBS  (the  only  available  experimental  data  set  for  this  case)  are  compared  with  the  UM  RANS 
computation  and  with  the  AFRL  VPM  2D  flat  plate  computation,  UM  RANS  2D  flat  plate  computation 
and  2D  quasi-steady  result.  The  latter  three  are  all  divided  by  2,  in  the  spirit  of  quasi-steady  thin  airfoil  or 
slender  body  theories,  which  both  predict  a  lift  curve  slope  of  na.  Rescaling  of  the  UM  2D  plate  RANS 
computation  gives  a  very  different  trend  from  the  AR  =  2  UM  RANS  computation.  Rescaling  of  the  AFRL 
VPM  2D  plate  computation  under-predicts  the  lift  coefficient  variation  and  produces  a  phase  lag  missing 
in  the  other  results,  which  appears  as  a  hysteresis  in  the  plot  with  respect  to  angle  of  attack.  The  other 
three  curves  -  AR  =  2  UM  RANS,  the  rescaled  quasi-steady  solution,  and  the  TUBS  measurement  -  are 
all  very  close,  showing  very  small  (if  any)  phase  lag  and  therefore  very  minor  hysteresis.  The  resulting  lift 
coefficient  time  history  therefore  has  the  appearance  of  behaviour  fully  consistent  with  linear  theory: 
quasi-steady  and  with  lift  curve  slope  according  to  slender-body  or  lifting-line  theory. 


Figure  6-9:  Lift  Coefficient  for  AR  =  2  Pitching-Piunging  Piate; 
Time  History  (ieft)  and  vs.  Angie  of  Attack  (right). 


Drag  coefficient  results  for  the  pitching-plunging  case  are  shown  in  Figure  6-10.  The  TUBS  measurement 
has  what  looks  like  a  steady  offset  of  Cd  ~  0.05  above  the  two  available  computations,  both  from  UM. 
Unlike  for  lift,  the  computed  drag  coefficient  for  the  wall-to-wall  plate  is  not  divided  by  2  to  overlay  with 
the  AR  =  2  result,  but  is  taken  directly  as  the  infinite-span  value.  Of  course,  the  entire  scaling  argument  for 
drag  is  rather  specious,  as  the  AR  =  2  plate  will  have  induced  drag  that  is  not  present  in  2D.  We  can  not, 
at  this  point,  make  a  definitive  argument  in  trends  in  drag  coefficient.  There  is  however  a  thrust  spike  near 
t/T  ~  0.25,  evidently  from  the  LEV.  The  mean  streamwise  force  is  a  net  drag. 
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Figure  6-10:  Drag  Coefficient  for  AR  =  2  Pitching-Piunging  Piate; 
Time  History  (ieft)  and  vs.  Angie  of  Attack  (right). 


6.4.2  Forces  for  the  AR  =  2  Plate  in  Pure-Plunge 

Lift  and  drag  coefficients  for  the  pure-plunging  case  are  given  in  Figure  6-11  and  Figure  6-12,  respectively. 
For  pure-plunge,  two  experimental  data  sets  are  available:  TUBS  and  UF  REEF,  both  from  wind  tunnels. 
The  intuitive  desire  is  to  remark  that  the  trends  in  the  SD7003  airfoil  and  wall-to-wall  cases  are  followed  for 
the  AR  =  2  plate,  in  that  force  coefficient  variations  for  the  pure-plunge  kinematics  are  less  than  for  pitch- 
plunge,  when  scaled  by  the  peak-to-peak  coefficient  variation.  While  this  assertion  is  not  really  false  for  the 
AR  =  2  plate,  it  is  largely  moot,  because  both  pitch-plunge  and  pure-plunge  kinematics  show  such  close 
agreement  amongst  the  various  data  sets.  As  with  the  pitch-plunge  case,  the  outliers  are  VPM,  and  in  the 
bottom  half  of  the  downstroke,  the  Li-scaling  of  the  UM  2D  RANS.  The  scaled  UM  2D  RANS  actually 
performs  better  for  pure-plunge  than  pitch-plunge,  so  in  that  sense  the  above-mentioned  assertion  remains 
relevant.  All  of  the  other  results  are  mutually  very  close  with  little  or  no  phase  lag,  and  follow  the  na  lift 
curve  slope  closely. 
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Figure  6-1 1 :  Lift  Coefficient  for  AR  =  2  Piate  in  Pure-Piunge; 
Time  History  (ieft)  and  vs.  Angie  of  Attack  (right). 


Figure  6-12:  Drag  Coefficient  for  AR  =  2  Piate  in  Pure-Piunge; 
Time  History  (ieft)  and  vs.  Angie  of  Attack  (right). 


In  drag  coefficient,  the  TUBS  measurement  has  a  negative  (thrusting)  offset  relative  to  the  REEF 
measurement  and  UM  AR  =  2  RANS;  as  in  the  pitch-plunge  case,  the  TUBS  drag  offset  is  on  the  order  of 
0.05,  but  in  the  opposite  direction.  Hysteresis  is  again  low,  and  the  drag  curve  is  much  more  sinusoidal  for 
pure-plunge  than  for  pitch-plunge.  Indeed,  we  find  that  AR  =  2  plate  results  are  more  quasi-steady  than  for 
the  wall-to  wall  (or  2D)  plate.  The  rescaled  UM  wall-to-wall  flat  plate  RANS  tracks  the  AR  =  2  result  for 
t/T  >  ~  0.45,  where  the  flow  is  largely  attached,  but  grossly  overshoots  earlier  in  the  plunge  downstroke, 
where  evidently  there  is  a  marked  difference  between  2D  and  AR  =  2  LEV  evolution. 
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6.5  EXPERIMENTS  ON  AR  =  2  PLATE  PURE-PLUNGE  AT  SMALLER 
STROUHAL  NUMBER 

As  with  the  SD7003  airfoil  and  wall-to-wall  plate,  it  remains  to  run  extended  parameter  studies  on  reduced 
frequency,  reduced  amplitude  and  so  forth.  One  of  our  principal  themes  is  the  extent  of  validity  of  quasi¬ 
steady,  linear  assumptions  to  unsteady  problems  with  flow  separation.  If  there  is  an  effective  angle  of 
attack  and  force  coefficients  depend  on  that  angle  of  attack,  then  the  actual  kinematics  of  how  that  angle  of 
attack  was  achieved  are  unimportant.  Otherwise,  if  the  kinematics  inform  the  force  coefficient  time  history 
directly,  then  the  very  concept  of  angle  of  attack  breaks  down,  and  with  it  superposition  and  so  forth. 
One  route  of  investigation,  already  suggested,  is  to  vary  the  motion  kinematics  such  that  the  reduced 
frequency  k  and  effective  angle  of  attack  history  (in  the  quasi-steady  sense)  are  kept  constant,  by  varying 
h  and  k.  This  in  fact  was  the  approach  taken  by  the  DLR  experimental  group.  The  idea  here  is  to  consider  a 
pure-plunge  motion  {X  =  0),  but  with  effective  angle  of  attack  history  the  same  or  nearly  the  same  as  that 
of  the  pitch-plunge  motion  {X  =  0.6).  For  k  =  0.25  this  results  mh  =  0.206. 

PIV  streamwise-component  velocity  contour  plots  are  given  in  Figure  6-13  for  four  phases  of  motion 
(top  of  plunge  stroke,  halfway  down,  bottom  of  stroke  and  halfway  up).  Regions  in  the  flowfield  plagued 
by  laser  reflections  from  the  model  surface  are  blanked-out. 
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Figure  6-13:  DLR  PIV;  /c  =  0.25  h  =  0.206  Pure-Plunge  of  AR  =  2  Plate,  with  a(t)  -  Pitch-Plunge 
of  .A  =  0.6;  Phases  Phi  =  0  (top  left),  90  (top  right),  180  (bottom  left)  and  270  (bottom  right). 


An  alternative  to  the  above,  still  retaining  the  same  nominal  angle  of  attack  time  history,  is  to  reduce  ^  by  a 
factor  of  2,  while  doubling  h  (Figure  6-14).  Strouhal  number  is  thus  also  unchanged.  Streamwise  velocity 
contours  for  this  case  are  in.  Qualitative  similarity  between  the  two  cases  is  excellent,  but  unfortunately 
detailed  comparison  is  precluded  by  reflections  from  the  model  surface.  Direct  force  measurements  are  not 
available.  Nevertheless,  we  have  good  supporting  evidence  in  all  of  our  cases  that  even  moderate  agreement 
between  flowfields  implies  excellent  agreement  in  lift  time  history  (but  certainly  not  vice  versa!).  We  can 
conclude,  at  least  preliminarily,  that  superposition  is  tenably  valid  in  the  sense  of  retaining  a  concept  of 
effective  angle  of  attack  through  various  combinations  of  reduced  frequency,  reduced  amplitude  and  X. 
This,  in  turn,  should  have  favourable  implications  for  designers  of  MAV  flapping  wings,  at  least  for  motions 
in  the  frequency  and  angle  of  attack  range  considered  so  far. 
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Figure  6-14:  DLR  PIV;  /c  =  0.125  A?  =  0.412  Pure-Plunge  of  AR  =  2  Plate,  with  a(t)  --  Pitch-Plunge 
of  A  =  0.6;  Phases  Phi  =  0  (top  left),  90  (top  right),  180  (bottom  left)  and  270  (bottom  right). 

Next  we  turn  to  an  entirely  different  problem,  returning  to  the  SD7003  airfoil,  with  an  even  lower 
amplitude  but  a  much  higher  frequency.  Here  too  we  wish  to  investigate  validity  of  quasi-steady  concepts 
and  traditional  small-disturbance  approximations  in  aerodynamics. 
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Chapter  7  -  HIGH-FREQUENCY  LOW-AMPLITUDE 
AIRFOIL  PURE-PLUNGE  OSCILLATIONS 

7.1  MOTIVATIONS 

So  far  we  have  considered  but  a  single  reduced  frequency,  k  =  0.25.  This  is  likely  representative  of  the 
flapping-flight  of  birds  in  cruise,  and  thus  presumably  of  flapping-wing  MAVs  of  the  larger  scale. 
And  this  is  towards  the  extreme  upper  bound  of  reduced  frequencies  associated  with  classical  dynamic 
stall.  It  is  doubtless  worthwhile  to  undertake  a  broad  parameter  study,  say  from  k  =  0.05  through  0.25,  1.0 
and  higher,  but  of  course  practical  limitations  such  as  lack  of  fanatically  studious  graduate  students 
preclude  this.  As  a  partial  effort,  we  consider  a  very  high  reduced  frequency,  relevant  perhaps  to  the 
flapping-flight  of  smaller  MAVs  operating  at  speeds  just  above  hover. 

The  k  =  3.93,  h  =  0.05  (St  =  0.125)  pure-plunge  was  considered  in  a  detailed  study  by  Jones  et  al.  [21]  on 
the  NACA  0012,  at  Re  -  10000,  as  part  of  a  parameter  study  of  thrust  efficiency  and  wake  topology  vs. 
Strouhal  number.  We  consider  this  problem,  also  in  pure-plunge,  but  for  the  SD7003,  at  Re  =  40000  and  at 
Re  =  10000,  and  an  offset  angle  of  attack  of  4°.  Again  we  compare  experiments  and  computations. 
As  the  motion  is  now  too  aggressive  for  wind  tunnels,  only  a  water  tunnel  run  is  considered  on  the 
experimental  side,  while  for  computation  we  again  compare  a  vortex  particle  method,  RANS  and  LES. 

7.2  FLOWFIELD  RESULTS 

In  the  following,  four  phases  of  the  plunge  sinusoidal  motion  are  reported.  Phase  phi  =  0  is  halfway  down 
on  the  plunge  stroke,  where  the  effective  angle  of  attack  is  largest  positive;  phi  =  90  is  at  the  bottom  of  the 
downstroke,  where  plunge-induced  angle  of  attack  is  zero;  phi  =  180  is  halfway  on  the  upstroke,  where  the 
effective  angle  of  attack  is  maximally  negative;  and  finally,  phi  =  270  is  the  top  of  the  plunge  stroke, 
where  again  plunge-induced  angle  of  attack  is  zero. 

AFRL  water  tunnel  instantaneous  and  phase-averaged  PIV  measurements  for  vorticity  are  presented  in 
Figure  7-1,  and  Figure  7-2  gives  phase-averaged  streamwise  component  of  velocity.  Phase-averages  are 
based  on  115  PIV  image  pairs,  from  a  sequence  of  120;  the  first  5  pairs  are  neglected  to  avoid  start-up 
transients  [97]. 
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Figure  7-1:  AFRL  PIV;  Re  =  40  K;  Vorticity  Contours,  k=  3.93  Plunge;  Phases  Phi  =  0, 
90, 180  and  270;  Phase-Averaged  Measurements  (left  column)  and 
Representative  Instantaneous  Measurements  (right  column). 


Figure  7-2:  AFRL  PIV;  Re  =  40  K;  Velocity  Contours,  k=  3.93  Plunge;  Phases 
Phi  =  0  (top  left),  90  (top  right),  180  (bottom  left)  and  270  (bottom  right). 


LES  computations  for  instantaneous  and  phase-averaged  vorticity,  and  for  phase-averaged  velocity,  are  in 
Figure  7-3  and  Figure  7-4,  respectively.  Instantaneous  and  phase-averaged  results  show  the  same  general 
flowfield  structure  away  from  the  airfoil  boundary  layer,  though  of  course  the  details  of  turbulence  will  be 
smeared  in  the  average.  Indeed  the  flow  is  transitional,  in  the  sense  that  transition  is  associable  with  the 
wall-bounded  small  LEV  visible  in  Figure  7-1  and  Figure  7-3.  The  Re  =  10  K  case  was  also  considered; 
details  are  reported  by  Visbal  [98],  and  for  brevity  the  present  discussion  is  limited  to  instantaneous  vs. 
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averaged  vorticity  contours  at  phi  =  0  (Figure  7-5).  There  transition  does  not  occur  until  the  trailing  edge 
region  and  the  near-wake,  and  the  wall-bound  LEV  is  a  discrete,  homogenous  blob  of  vorticity  that  differs 
little  between  instantaneous  and  phase-averaged.  The  qualitative  conclusion  is  that  this  flow  is  highly 
periodic  at  both  Reynolds  numbers,  with  a  reverse  Karman  vortex  street,  and  is  therefore  much  akin  to  the 
findings  for  the  NACA0012  at  zero  mean  angle  of  attack  [21].  In  turn,  this  implies  insensitivity  to  the 
airfoil  shape  and  the  Reynolds  number. 


Figure  7-3:  AFRL  LES  Computations;  Re  =  40  K;  Vorticity  Contours,  k=  3.93  Piunge;  Phases 
Phi  =  0,  90,  180  and  270;  Phase-Averaged  Measurements  (ieft  coiumn)  and 
Representative  instantaneous  Measurements  (right  coiumn). 
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Figure  7-4:  AFRL  LES  Computations;  Re  =  40  K,  Veiocity  Contours,  k-  3.93  Piunge; 
Phases  Phi  =  0  (top  ieft),  90  (top  right),  180  (bottom  ieft)  and  270  (bottom  right). 


Figure  7-5:  AFRL  LES  Computations;  Re  =  10  K;  Vorticity  Contours,  k-  3.93  Piunge; 
Phase  Phi  =  0;  Phase-Averaged  (top  ieft)  and  instantaneous  (top  right); 
and  Phase-Averaged  AFRL  PiV  (bottom). 


Next  we  turn  to  RANS  at  Re  =  40  K  (University  of  Michigan).  Velocity  contours  are  shown  in  Figure  7-6, 
and  vorticity  is  in  Figure  7-7.  Agreement  with  phase-averaged  LES  is  very  close,  which  is  to  say  that 
RANS  without  a  transition  model  is  adequate  for  capturing  the  overall  wake  structure  for  this  unsteady 
transitional  periodic  flow,  in  the  sense  of  a  phase-average. 
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Figure  7-6:  UM  RANS  Computations;  Re  =  40  K;  Veiocity  Contours,  k-  3.93  Piunge; 
Phases  Phi  =  0  (top  ieft),  90  (top  right),  180  (bottom  ieft)  and  270  (bottom  right). 


Figure  7-7:  UM  RANS  Computations;  Re  =  40  K;  Vorticity  Contours,  k-  3.93  Piunge; 
Phases  Phi  =  0  (top  ieft),  90  (top  right),  180  (bottom  ieft)  and  270  (bottom  right). 


Finally,  we  compare  with  a  2D  vortex-particle  calculation  (AFRL  VPM).  Traces  of  vortex  particles  released 
at  the  airfoil  TE  are  shown  in  Figure  7-8,  which  are  analogous  to  PIV  vorticity  in  Figure  7-1  and  LES 
vorticity  in  Figure  7-3;  roughly,  spatial  concentration  of  vortex  particles  is  akin  to  local  intensity  of  vorticity. 
The  present  implementation  of  the  VPM  algorithm  has  no  provision  for  shedding  vorticity  from  the  airfoil 
anywhere  except  for  the  trailing  edge,  in  particular  missing  the  leading-edge  rollers  on  the  airfoil  suction 
surface.  Nevertheless,  the  wake  time  history  is  very  similar  to  those  of  the  LES  computation,  the  experiment 
and  the  RANS.  It  appears,  therefore,  that  LEV  production  for  this  kind  of  flow  is  unimportant  for  the  overall 
wake  structure. 
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Figure  7-8:  Vortex  Particle  Calculation  for  k  =  3.93  Plunge;  Phases  Phi  =  0  (top  left), 
90  (top  right),  180  (bottom  left)  and  270  (bottom  right). 


In  sum,  then,  all  of  the  computational  methods  give  similar  wake  structure,  though  the  fidelity  of  resolution 
of  boundary  layer  structure  goes  in  proportion  to  the  complexity  of  the  method.  Following  the  by  now 
established  theme  of  our  work,  we  next  consider  how  the  calculation  of  aerodynamic  force  compares 
amongst  the  various  methods;  if  the  computed  flowfields  are  similar,  the  computed  loads,  if  past 
performance  is  any  guarantee  of  future  results,  should  be  even  more  similar. 


7.3  FORCE  COMPUTATIONS 

Lift  coefficient  time  history  from  LES  computation  at  Re  =  40  K  and  10  K,  RANS  computations  at 
Re  =  40  K  and  from  Theodorsen’s  formula,  is  given  in  Figure  7-9.  Reynolds  number  seems  to  have 
essentially  no  effect  between  the  two  LES  runs,  despite  the  difference  in  boundary  layer  properties  and 
role  of  transition.  RANS  splits  the  difference  between  LES  and  Theodorsen,  which  itself  is  not  large. 
Thus  Theodorsen’s  formula  gives  acceptable  lift  time  history  prediction  despite  the  wake  being  far  from 
planar,  and  the  obvious  presence  of  vortex  shedding.  This  seems  counterintuitive;  why  would  Theodorsen’s 
formula  work  despite  the  blatant  violation  of  its  underlying  assumptions?  The  resolution  is  in  the  fact  that 
non-circulatory  force  dominates  circulatory  force,  as  seen  from  comparing  just  the  non-circulatory  portion  of 
Theodorsen’s  formula  with  the  LES  solution  (green  curve  in  Figure  7-9).  Theodorsen’s  equation  is  valid  for 
arbitrarily  large  frequency,  but  only  small  amplitude.  Evidently,  despite  the  reverse  Karman  vortex  street, 
h  =  0.05  is  a  “small  enough”  amplitude.  Linear  superposition  also  seems  to  be  valid  here. 
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AFRL  LK^  k  3  113  Re  40k 


Figure  7-9:  Lift  Coefficient  Time  History  (ieft)  and  vs.  Angie  of  Attack  (right); 
k  =  3.93  h  =  0.05  Pure-Piunge  of  SD7003  Airfoii  at  Re  =  40  K  and  10  K. 


Drag  coefficient  is  shown  in  Figure  7-10,  and  integrated  lift  and  drag  over  the  plunge  cycle  are  tabulated 
in  Table  7-1.  Consistent  with  the  reverse  Karman  vortex  wake,  a  net  thrust  (negative  drag)  is  produced. 
However,  as  Re  decreases,  evidently  viscous  drag  increases,  and  the  net  thrust  is  attenuated. 


Figure  7-10:  Drag  Coefficient  Time  History;  k=  3.93  h  =  0.05 
Pure-Piunge  of  SD7003  Airfoii  at  Re  =  40  K  and  10  K. 
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Table  7-1:  Computed  Integrated  Force  Coefficients;  k  =  3.93  h  =  0.05  Pure-Plunge 


Cl 

Cd 

LES  40  K 

0.73 

-0.081 

LES  10  K 

0.66 

-0.054 

RANS  40  K 

0.58 

-0.069 

It  remains  to  explore  higher  amplitudes  to  see  where  Theodorsen’s  formula  convincingly  breaks  down, 
and  for  that  matter  to  find  clear  distinction  between  the  results  of  vortex  particle  methods,  RANS  and  LES. 
It  also  remains  to  consider  start-up  transients  at  the  commencement  of  a  periodic  motion,  or  alternatively 
non-periodic  problems  such  as  pitch  ramp-and-hold  or  ramp-and-return. 

We  conclude,  at  least  tentatively,  that  high-frequency  low-amplitude  cases  are  “easy”,  in  the  sense  that  the 
imposed  motion  dynamics  constrains  the  aerodynamic  solution  to  a  periodicity  amenable  even  to  inviscid 
analysis,  and  at  high  enough  reduced  frequency  the  lift  coefficient  is  dominated  by  the  so-called  non- 
circulatory  loads.  Of  course,  it  is  the  fluid  streamlines  that  determine  the  pressure  gradients  responsible  for 
non-circulatory  loads,  and  the  shed  and  bound  vortex  strength,  together  with  the  trailing  edge  Kutta  condition 
and  no-flow-through  boundary  condition  determine  the  streamlines.  Thus  all  aerodynamic  forces  are  in  the 
end  “circulatory”.  The  point,  however,  is  that  one  may  blithely  disregard  shed  and  even  bound  vortex 
strength,  retaining  in  the  equations  only  the  terms  related  to  fluid  acceleration,  and  still  obtain  a  workable 
approximation  -  if  the  reduced  frequency  is  high  enough.  We  leave  as  a  future  problem  to  ascertain  how  far 
this  can  be  extended  to  moderate  or  large  motion  amplitudes,  and  to  non-periodic  motions. 

Thus  far  our  work  has  been  exclusively  on  rectilinear  or  longitudinal-plane  motions,  where  the  airfoil  or 
plate  oscillates  in  a  plane  parallel  to  the  free-stream  and  normal  to  the  plane  of  the  wing  planform.  This  is 
the  most  straightforward  abstraction  for  the  flight  of  large  birds,  for  example  -  but  not  necessarily  for  the 
smaller  flying  creatures  such  as  insects  or  hummingbirds.  While  a  collective  treatment  across  multiple 
laboratories  has  in  the  present  work  been  intractable  for  non-rectilinear  motions,  we  present  in  the  next 
chapter  a  promising  line  of  investigation  for  a  rotating  or  “whirling”  motion. 
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Chapter  8  -  NON-RECTILINEAR  MOTIONS: 

WHIRLING  FLAT  PLATE 

8.1  PURPOSE  AND  TYPE  OF  ^WHIRLING”  MOTIONS 

8.1.1  Rotational  Motions  and  Leading  Edge  Vortex  Retention 

There  is  presently  considerable  controversy  about  the  stability  of  insect-flapping  LEVs.  As  discussed  above, 
many  mechanisms  for  LEV  stability  have  been  proposed,  and  indeed  the  whole  issue  of  LEV  stability  is  but 
one  of  several  mechanisms  for  explaining  the  unexpectedly  high  lift  coefficients  encountered  in  insect 
flapping-wing  studies.  One  of  the  main  controversies  is  whether  there  is  something  intrinsic  to  “flapping”  in 
the  sense  of  pivoting  or  swinging  about  a  fixed-point  at  the  wing  root,  as  opposed  to  rectilinear  translation, 
which  was  the  subject  of  all  of  the  preceding  chapters  in  this  report. 

A  quasi-steady  setup  developed  by  Usherwood  and  Ellington  rotates  a  wing  like  a  propeller  in  a  water  tank 
with  no  ambient  flow  [99].  When  a  wing  undergoes  a  whirling  propeller-like  motion  there  is  a  velocity 
gradient  from  the  wing  root  to  tip  not  present  in  rectilinear  motion.  The  propeller  experiments  preserve  the 
three-dimensional  flow  as  a  velocity  gradient,  but  neglect  the  starting  and  stopping  that  occurs  at  the 
beginning  and  end  of  the  wing  stroke.  The  idea  is  to  observe  the  long-term  properties  of  leading  edge  vortex 
formation,  and  retention  or  shedding;  once  formed,  does  the  LEV  stay  resident  near  the  wing  leading  edge 
for  long  dimensionless  time,  or  does  it  detach  and  convect  away  from  the  wing,  as  we  saw  in  the  rectilinear 
periodic  oscillations  in  previous  chapters?  In  brief,  the  conclusion  of  the  propeller-like  experiments  is  that 
long-term  LEV  retention  is  possible,  at  least  at  the  Reynolds  numbers  of  the  subject  experiment  [99]. 

But  as  an  insect  wing  stroke  is  only  2  to  5  chord  lengths  [100,101],  unsteady  effects  at  the  beginning  and 
end  of  the  flapping  cycle  can  produce  significant  additional  lift  and  raise  questions  beyond  the  propeller- 
type  experiments  on  LEV  formation/retention  history  and  transient  effects.  One  is  interested  in  the  role  of 
starting  transients  -  how  acceleration  profile  and  intensity  affects  the  subsequent  fiowfield  and  lift 
coefficient  evolution,  and  how  long  it  takes  to  relax  to  quasi-steady  (if  one  exists)  state.  This,  of  course, 
was  ignored  in  our  work  reported  above,  where  the  experiments  and  computations  were  deliberately  run 
long  enough  for  starting  transients  to  relax,  so  that  phase-averages  would  be  taken  over  at  least  nominally 
periodic  conditions. 

8.1.2  Details  of  Motion  Kinematics 

The  waving  wing  experiment,  described  in  an  earlier  chapter,  is  designed  to  mimic  the  translational  phase 
of  an  insect  wing  stroke,  building  on  Ellington-type  propeller  experiments  to  include  wing  starting  and 
stopping.  Three  angular  velocity  starting  profiles  were  considered:  linear,  sinusoidal,  and  exponential, 
each  reaching  the  maximum  velocity  necessary  to  achieve  a  Reynolds  number  of  60000  in  0.10,  0.25, 
or  0.60  chord-lengths  of  travel  referenced  at  3/4  span.  The  commanded  velocity  and  acceleration  profiles 
are  shown  in  Figure  8-1  as  plotted  vs.  time,  and  in  Figure  8-2  as  plotted  vs.  angular  displacement. 
An  exponential  rise  in  rotation  rate  vs.  time  produces  a  linear  rise  in  rotation  rate  vs.  angular  displacement. 
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Figure  8-1:  Commanded  Wing  Kinematics  Piotted  vs.  Time,  for  First  Quarter-Stroke  Acceierating 
Over  0.25  c;  Anguiar  Veiocity  (soiid)  and  Acceieration  (dashed)  as  a  Function  of  Time. 


Figure  8-2:  Commanded  Wing  Kinematics  Piotted  vs.  Anguiar  Dispiacement,  for  the  First 
Quarter-Stroke;  Anguiar  Veiocity  for  a  Wing  Acceierating  Over  the  First  0.25  c  of  Travei. 


8.2  FORCE  DATA 

Figure  8-3  shows  the  averaged  and  smoothed  lift  coef  f  icients  for  a  waving  wing  at  a  5  and  15  degree  angle 
of  attack  accelerating  to  maximum  velocity  in  0.10,  0.25  and  0.60  chords.  It  can  be  seen  that  changing  the 
angle  of  attack  primarily  changes  the  magnitude  of  the  measured  force  while  the  shape  of  the  curve  remains 
much  the  same.  All  of  the  curves  are  characterized  by  an  initial  steep  increase  in  lift  coefficient  leading  to  a 
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lift  peak,  followed  by  a  sharp  drop  and  subsequent  recovery  to  an  intermediate  level  that  is  then  maintained 
for  the  rest  of  the  constant  velocity  portion  of  the  wing  stroke.  It  appears  that  after  2  to  3  chord-lengths  of 
travel  (at  3/4  span,  corresponding  to  34  to  51  degrees  of  rotation)  the  lift  nears  a  steady-state  value. 
The  shape  of  the  lift  curve  shown  here,  including  the  initial  lift  peak,  is  very  similar  to  those  reported  by 
Dickinson  and  Gotz  at  Re  <  1,000  [102]. 
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Figure  8-3:  Lift  Coefficients  for  a  Waving  Wing  Acceierating  Over  0.10  c  (top  ieft),  0.25  c 
(top  right),  0.60  c  (bottom  ieft)  and  Aii  Exponentiai  Veiocity  Profiies  (bottom  right). 


The  main  differences  due  to  different  acceleration  patterns  (Figure  8-4)  are  observed  in  the  early  stages 
of  the  wing  stroke.  The  slowest  acceleration  profile  is  the  exponential  over  0.60  chords,  not  even  visible  in 
Figure  8-4,  with  the  next  slowest  being  the  constant  acceleration.  Both  of  these  velocity  profiles  result  in  a 
relatively  low  delayed  lift  peak  past  1  chord  of  travel  as  in  Figure  8-3.  Maximum  lift  coefficients  for  the 
linear,  sinusoidal,  and  exponential  velocity  profiles  accelerating  in  0.60  c  are  0.87,  0.88  and  1.01. 
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Figure  8-4:  Angular  Acceleration  for  the  Three  Velocity  Profiles:  Linear  (solid  line),  Sinusoidal  (dashed  line) 
and  Exponential  (dash-dotted  line)  Over  0.10  c  (left),  0.25  c  (center)  and  0.60  c  (right). 


At  higher  accelerations  the  effect  of  acceleration  on  the  lift  peak  diminishes  and  all  of  the  curves  show  a 
pronounced  maximum  approximately  50%  above  the  steady  state  value  near  0.60  chord-lengths  travelled, 
well  into  the  constant  velocity  portion  of  the  wing  stroke.  The  largest  acceleration  studied  is  the  exponential 
over  0.10  chords,  with  a  maximum  of  40.2  rad/s  and  a  maximum  lift  coefficient  of  1.42.  Over  0.25  chords 
the  maximum  acceleration  is  8.6,  13.6,  and  1.6  rad/s  for  the  linear,  sinusoidal,  and  exponential  velocity 
profiles  with  maximum  lift  coefficients  of  1.34,  1.39,  and  1.26,  respectively,  as  in  Figure  8-2.  Despite  the 
large  difference  in  the  acceleration,  the  lift  is  of  similar  magnitude  and  timing.  It  is  interesting  to  note  that  the 
timing  of  the  transient  lift  peak  does  not  change  significantly  whether  the  wing  is  accelerated  over  0.10  or 
0.25  chords,  suggesting  that  the  sharp  lift  peak  evident  in  the  force  data  is  due  to  aerodynamic  effects  rather 
than  inertia.  Lift  curves  at  a  5  degree  angle  of  attack  look  much  like  those  at  15  deg  with  a  lower  CLmax  and 
“steady-state”  value.  Cl  levels  off  near  0.95  for  a  =  15  deg  or  0.27  for  a  =  5  deg  after  2.5  chords  of  travel  for 
all  wing  kinematics.  For  any  profile  accelerating  in  0.10  or  0.25  chords,  the  maximum  Cl  occurs  near 
0.6  chords  of  travel.  Beckwith  and  Babinsky  [103]  performed  experiments  on  an  impulsively  started  fiat 
plate  wing  towed  in  pure  translation  through  a  water  tank.  Figure  8-5  shows  the  unsteady  lift  coefficients 
measured  on  the  waving  wing  as  compared  to  their  sliding  wing.  Only  the  constant  acceleration  over  0.60  c 
is  shown  here  as  faster  accelerations  are  not  available  for  the  sliding  wing.  The  force  history  for  the  two 
cases  is  very  similar.  Also  shown  is  the  expected  lift  according  to  Wagner’s  theory  using  the  long-term 
steady  state  value  measured  on  the  translating  wing  as  the  asymptotic  lift  coefficient.  At  a  5  degree  angle  of 
attack,  below  steady  state  stall,  the  unsteady  development  of  lift  coefficient  for  both  experiments  closely 
matches  those  predicted  by  Wagner  [28].  At  larger  angles  of  attack  and  in  the  initial  phase  of  the  wing  stroke 
other  aerodynamic  effects  appear  to  dominate  the  force  history.  Lift  curves  at  the  higher  angle  of  attack 
resemble  those  at  Reynolds  numbers  less  than  1,000  reported  by  Dickinson  and  Gotz  [102]. 
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Figure  8-5:  Comparison  of  Unsteady  Lift  Coefficients  Observed  for 
a  Waving  Wing  and  Transiating  Wing  Acceierating  Over  0.60  c. 

Cl  versus  a  for  the  waving  wing  (constant  acceleration  over  0.25  c)  at  an  angle  of  attack  from  0  to  25 
degrees  are  given  in  Figure  8-6.  Here,  “max  Cl”  is  the  maximum  Cl  achieved  at  the  lift  peak  as  shown  in 
Figure  8-2.  The  steady  Cl  refers  to  the  mean  value  measured  between  3  and  4  chords  of  travel  where  lift 
levels  out.  Both  lift  coefficients  appear  to  be  linear  in  a.  The  line  of  best  fit  for  the  maximum  lift  coefficient 
is  1.82  na  and  1.31  na  for  the  quasi-steady  lift  coefficient.  For  comparison,  Wagner  predicts  1.6  na  after 
3.5  chords  of  travel.  Also  shown  is  the  2  na  curve  from  thin  airfoil  theory  and  Theodorsen’s  theory  [26] 
where  the  Theodorsen  function,  C(k),  is  C(k)  =  0.52  -  0.06  i  for  acceleration  over  0.25  c.  The  maximum  lift 
achieved  on  the  waving  wing  is  similar  to  the  thin  airfoil  theory  prediction  despite  the  fundamental 
differences  in  the  flow-field,  especially  at  higher  angles  of  attack.  The  extension  of  Theodorsen’s  theory  to 
the  waving  wing  under-predicts  lift  at  all  of  the  tested  angles  of  attack.  Beckwith  and  Babinsky’s 
experiments  on  a  sliding  wing  suggest  that  the  true  steady  state  lift  is  not  achieved  until  the  wing  has 
travelled  much  further  and  is  lower  than  the  lift  after  3  to  4  chord-lengths  of  travel  [103].  This  may  explain 
the  under-prediction  shown  here. 
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a,  deg 


Figure  8-6:  Lift  Coefficient  as  a  Function  of  Angie  of  Attack. 


8.3  FLOWFIELD  MEASUREMENTS  WITH  PARTICLE  IMAGE 
VELOCIMETRY 

8.3.1  Vortex  Tracking 

Considering  the  flow  during  the  early  stages  of  motion,  Figure  8-7  compares  both  the  velocity  field  and  A 
contours  to  the  force  history  at  three  points  in  time  for  the  waving  wing  undergoing  constant  acceleration 
over  0.25  chord-lengths  of  travel  at  %  span.  is  a  scalar  proxy  for  vorticity,  which  is  useful  in  complex 
flows  to  isolate  rotational  effects  from  small  concentrations  of  shearing  and  other  perturbations  [104], 
and  contours  of  A  are  thus  a  form  of  vortex  tracking.  There  are  noticeable  differences  in  the  flow  pattern  at 
these  points  in  time.  First,  (a)  a  distinct  leading  edge  vortex  forms  on  the  suction  side  of  the  plate. 
The  starting  vortex  is  visible  just  downstream  of  the  trailing  edge.  At  the  minimum  lift  point  (b), 
this  vortex  has  shed  and  moved  downstream  while  a  new  leading  edge  vortex  forms.  As  the  lift  approaches 
the  steady-state  value  (c)  several  vortices  coexist  above  the  wing.  After  this  time  there  is  a  continuous 
pattern  of  shedding  vortices  in  a  quasi-periodic  fashion.  The  images  in  Figure  8-7  suggest  that  the  rapid 
build-up  of  lift  during  the  initial  stages  of  the  wing  stroke  is  due  to  the  development  of  the  attached 
leading  edge  vortex  and  the  movement  of  the  starting  vortex  away  from  the  wing. 
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Figure  8-7:  Flowfield  Observed  at  Three  Points  During  the  Wing  Stroke:  (a)  Near  the  Initiai  Peak 
at  x/c  =  0.45  c;  (b)  At  the  Lift  Minimum  at  x/c  =  1.50  c;  and  (c)  At  the  Onset  of  the  Quasi-Steady 
Period  x/c  =  2.39  c  -  Force  History  (left),  Velocity  Field  (centre)  and  Contours  of 
A  (right)  for  a  Whirling  Wing  with  Constant  Acceleration  Over  0.25  c. 

Recently  it  has  been  suggested  that  the  attachment  of  the  leading  edge  vortex  is  related  not  to  the  Reynolds 
number,  but  to  the  Rossby  number.  Lentink  and  Dickinson  observed  a  stable  LEV  at  low  Rossby  number 
(Ro  =  radius/chord)  over  a  range  of  Reynolds  numbers  for  a  propeller-like  motion  [105].  The  waving  wing 
experiment  has  been  performed  at  a  Rossby  number  of  both  2  and  4,  well  within  the  low  Ro  range 
identified  by  Lentink  and  Dickinson,  but  there  is  no  evidence  of  an  attached  leading  edge  vortex  over  an 
extended  period  of  time  on  the  waving  wing. 

8.3.2  Leading  Edge  Vortex  Circulation 

Again  considering  a  constant  acceleration  over  0.25  c,  the  circulation  of  the  LEV  was  computed  by 
integrating  along  a  contour  of  constant  A  to  find  the  circulation  in  the  usual  way,  E  =  ds. 

The  computed  circulation  was  then  normalized  by  the  local  freestream  velocity  and  wing  chord.  Figure  8-8 
shows  the  normalized  circulation  of  the  first  leading  edge  vortex  formed  at  half-span  at  a  25  deg  angle  of 
attack.  Each  data  point  on  the  chart  was  computed  from  a  separate  snapshot  of  the  velocity  field  as  measured 
using  PIV.  The  dashed  line  marks  x/c  =  0.58,  the  point  of  maximum  lift. 
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Figure  8-8:  Normalized  Circulation  Computed  from  PIV  Data  at  Half-Span;  AR  =  4,  a  =  25  deg. 


Figure  8-9  shows  contours  of  A>0.6  on  top  of  the  vorticity  field  for  an  aspect  ratio  4  wing  at  a  25  degree 
angle  of  attack  throughout  the  wing  stroke.  Each  column  represents  a  location  on  the  wing  span,  and  each 
row  represents  a  wing  stroke  angle.  From  the  start  of  the  wing  stroke,  an  attached  leading  edge  vortex 
begins  to  form  and  almost  all  of  the  vorticity  in  the  flowfield  is  contained  in  the  A  =  0.6  contour. 
The  circulation  of  this  leading  edge  vortex  grows  steadily  until  the  point  of  maximum  lift  at  x/c  =  0.58. 
From  Figure  8-8,  the  maximum  strength  of  the  attached  leading  edge  vortex  is  near  F*  =  1.3,  at  which 
point  the  leading  edge  vortex  has  elongated  and  a  second  region  of  high  vorticity  has  appeared  within  the 
A  =  0.6  contour.  As  the  wing  continues  to  move,  the  downstream  portion  of  the  leading  edge  vortex  sheds 
and  the  two  regions  of  vorticity  become  separate  vortex  cores. 
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(a)  b  ^  1/2,  e  ^  1 0.85  deg.  x/c  =  0.46 


(b)  b  ^  3/4,  e  ^  1 0.85  deg^  x/c  ^  0.65  (c)  b  ^  7/8,  8  =  10.85  deg.  x/c  ^  0.74 


(e)  b  =  3/4.  0  =  27.63  deg,  x/c  =  1 .65 


K 


(f)  b  =  7/8.  8  =  27.63  deg.  x/c  =  1.80 


(g)  b  -  I  /2. 9  -  40.00  deg.  x/c  -  1 .69 


(h)  b  -  3/4.  e  -  40.00  deg.  x/c  -  2.39 


(k)  b  -  3/4,  e  =  50.85  deg.,  x/c  =  3.04 


(i)  b  -  7/8. 0-40.00  deg,  x/c  -  2.74 


(I)  b  =  7/8.  0  =  50.85  deg,  x/c  =  3.48 


(n)  b  -  3/4.  e  -  67.63  deg.  x/c  -  4.05 


Figure  8-9:  Chordwise  View  of  Vorticity  (co/U'c)  and  A-Contours;  AR  =  4,  a  =  25  deg  - 
each  row  represents  the  three  stations  on  the  wing  at  an  instant  in  the  wing  stroke. 


Near  x/c  =  0.69  it  is  difficult  to  reliably  track  the  first  vortex  and  F*  values  jump  between  1.0  and  1.3. 
Most  of  the  vorticity  is  still  contained  within  the  higher-level  A-contours,  but  there  is  also  vorticity 
moving  between  the  two  vortex  cores.  As  the  wing  stroke  continues  the  separated  leading  edge  vortex 
continues  to  move  downstream  and  entrain  this  vorticity.  Eventually  this  vortex  rolls  up  at  into  a  more 
coherent  structure  at  x/c  =  0.89  and  the  F*  values  become  more  reliable.  The  first  vortex  continues  to 
move  downstream  and  away  from  the  still-attached  second  leading  edge  vortex.  Vorticity  is  still  present 
between  the  two  vortex  cores  and  the  first  vortex  continues  to  entrain  some  of  this  vorticity  causing  the 
vortex  to  strengthen  slightly. 

8.3.3  Three-Dimensional  Effects 

The  chordwise  view  presented  in  Figure  8-9  shows  the  vorticity  field  normalized  by  the  local  freestream 
velocity  in  order  to  reveal  variations  along  the  span.  It  can  be  seen  that  the  normalized  vorticity  is  higher  near 
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the  wing  root  (1/2  span)  than  further  outboard.  This  is  particularly  clear  in  (d)  and  (e).  Also,  the  leading  edge 
vortex  is  larger  near  the  wingtip  as  shown  in  (a)  -  (c). 

Figure  8-10  shows  the  spanwise  view  of  the  wing  as  the  trailing  edge  passes  through  the  light  sheet  at 
different  points  in  the  wing  stroke  (wing  waving  into  the  page,  suction  side  of  the  wing  on  the  right).  There  is 
a  clear  tip  vortex  visible  from  the  start  of  the  wing  stroke  (10.85  degrees)  and  a  shear  layer  that  becomes  very 
clear  at  67.63  degrees.  Late  in  the  wing  stroke,  the  tip  vortex  is  destroyed  by  the  vorticity  passing  over  the 
suction  side  of  the  wing.  At  lower  angles  of  attack  (5  and  15  degrees)  the  tip  vortex  persists  throughout  the 
wing  stroke. 


(a)  0  =  10  .S5  deg  (b)  0  =  27  .63  deg  (c)  9  =  40  .00  deg 


(d)  6  =  50  .85  deg  (e)  9  =  67  .63  deg  (0  9  =  80  .00  deg 

Figure  8-10:  Spanwise  View  of  Vorticity  (co/U'c)  and  A-Contours;  AR  =  4,  a  =  25  deg. 


8.3.4  Flow  Model  Development 

Based  on  these  observations,  we  propose  that  the  flow  on  an  impulsively  started  waving  wing  develops 
over  three  stages  as  illustrated  in  Figure  8-11.  At  the  start  of  the  wing  stroke,  during  the  initial  transient. 
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the  flow-field  is  characterized  by  the  growth  of  a  strong  leading  edge  vortex.  Flow  separates  at  the  sharp 
leading  edge  and  quickly  forms  a  leading  edge  vortex  that  remains  attached  to  the  wing.  During  this  phase 
the  strengthening  of  the  leading  edge  vortex  causes  a  rapid  increase  in  lift.  This  phase  ends  when  the 
vortex  sheds  from  the  leading  edge.  The  ultimate  strength  of  the  leading  edge  vortex  and  its  development 
is  influenced  by  the  wing’s  acceleration  pattern.  Lower  accelerations  produce  a  slower  vortex  growth  and 
reduce  the  strength  at  separation.  However,  above  a  certain  level  of  acceleration  the  vortex  development 
becomes  independent  from  the  kinematics. 


Estabiished  flow 


i\  \  2  ^  15 

irhiTiJ?  Irrrii.  li.'il 


Figure  8-11:  Proposed  Model  of  Three  Phases  of  Flowfield  Development  for  the  Whirling  Plate. 


The  next  flow  observed  on  the  waving  wing  is  the  developing  flow,  beginning  with  the  separation  of  the 
initial  leading  edge  vortex.  Flow  continues  to  separate  over  the  leading  edge  and  a  second  leading  edge 
vortex  forms.  Because  the  first  vortex  is  no  longer  close  to  the  wing  surface  its  effect  on  lift  is  diminished 
and  the  new  leading  edge  vortex  does  not  achieve  the  same  ultimate  strength  before  it  sheds.  In  this  phase 
the  total  lift  is  relatively  low. 

Finally,  in  the  established  flow  phase  a  periodic  pattern  of  vortex  shedding  from  the  leading  edge  is 
observed.  During  this  phase  there  are  typically  three  to  four  vortices  present  above  the  wing  at  all  times 
producing  a  relatively  high  lift  coefficient.  At  this  point  the  structure  of  the  flow-field  appears  to  settle  into 
a  periodic  vortex  shedding  but  it  is  not  yet  certain  how  many  chord-lengths  must  be  traveled  before  the 
flow  is  truly  in  a  periodic  quasi-steady  state. 


8.4  SUMMARY  OF  THE  WHIRLING  WING  EXPERIMENTS 

The  rotating  wing  experiment  is  a  fully  three-dimensional  simplification  of  the  flapping  wing  motion 
observed  in  nature  for  small  flyers  such  as  many  insects.  The  spanwise  velocity  gradient  and  stopping  and 
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starting  acceleration  preserve  key  features  of  the  wing  stroke  which  are  neglected  by  experiments 
concerned  with  established  periodic  flows. 

All  of  the  whirling- wing  flows  showed  similar  behavior.  High  lift  (approximately  1.5  times  the  quasi-steady 
value)  is  achieved  during  the  first  chord-length  of  travel  due  to  the  strong  attached  leading  edge  vortex.  After 
the  initial  leading  edge  vortex  separates,  lift  values  drop  as  the  LEV  passes  downstream  over  the  wing  and  a 
second  LEV  forms.  LEVs  continue  to  form,  grow,  and  separate,  but  subsequent  vortices  are  not  as  strong  as 
the  first  one  and  lift  values  level  out.  Because  the  translational  phase  of  the  typical  insect  wing  stroke  is 
between  2  and  4  chord-lengths  of  travel,  the  initial  lift  transient  and  subsequent  flow  development  phases  on 
a  waving  wing  are  likely  to  be  significant  for  lift  production  for  MAV  applications  inspired  by  insect-type  of 
flapping.  Wing  kinematics  had  only  a  small  effect  on  the  aerodynamic  forces  produced  by  the  rotating  wing. 
Quasi-steady  lift  and  drag  forces  were  very  similar  for  all  acceleration  patterns  tested.  In  the  early  stages  of 
the  wing  stroke  the  velocity  profile  affected  both  the  timing  and  magnitude  of  the  lift  peak  for  low 
accelerations.  However,  at  high  accelerations  the  velocity  profile  ceased  to  be  significant.  The  measured 
circulation  of  the  LEV  agrees  well  with  force  data  thus  supporting,  the  proposed  3 -stage  model  of  flow 
development. 

Of  course,  it  remains  to  study  whirling-wing  motions  at  the  high  angles  of  attack  more  typical  of  insect-type 
of  flight.  But  this  can  be  said  about  all  of  the  cases  in  our  work,  which  have  in  general  cleaved  towards 
classical-dynamic  stall  motivation.  As  we  conclude,  we  offer  recommendations  on  how  to  generalize 
investigations  both  in  extension  of  those  presented  above,  and  in  conceptually  new  directions. 
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Chapter  9  -  CONCLUSIONS  AND  OUTLOOK 

9.1  RECAPITULATION  OF  THE  TEST  CASES  IN  THE  PRESENT  STUDY 

There  is  a  vast  literature  in  unsteady  aerodynamics  at  low  Reynolds  number,  not  just  in  applications  to 
Micro  Air  Vehicles  or  natural  flyers,  but  for  airfoils  in  general,  for  aerodynamic  testing  of  aerospace 
configurations  in  small-scale  facilities  where  the  Reynolds  number  happens  to  be  low  not  by  intention  but 
by  necessity,  for  fundamental  studies  such  as  oscillating  circular  cylinders,  for  flow-control  experiments 
and  on  and  on.  Our  objectives  have  been: 

1)  To  extend  to  dynamic  problems  the  recent  work  on  static  low  Reynolds  number  airfoils,  laminar 
separation  bubbles,  experimental  and  computational  and  methods  for  such  airfoils; 

2)  To  benchmark  a  broad  range  of  experimental  facilities  and  computational  formulations  established 
for  research  programs  of  the  type  that  the  present  program  is  prototypical; 

3)  To  assess  the  maturity  and  utility  of  various  computational  methods  for  low-Re  unsteady 
aerodynamic  predictions; 

4)  To  assess  the  importance  of  flow  quality  in  wind  tunnels,  water  tunnels  and  tow  tanks; 

5)  To  observe  how  the  low-Re  unsteady  aerodynamics  problem  might  differ  from  the  classical  high- 
Re  problem  of  dynamic  stall;  and  most  importantly;  and 

6)  To  begin  to  bridge  the  gap  between  practical  MAV  designers  and  the  fluid  mechanics  community. 

Such  an  endeavour  is  necessarily  more  grandiose  in  planning  than  in  execution,  and  indeed  the  practical 
results  are  unavoidably  episodic.  We  are  limited  to  a  few  configurations  and  a  few  run  conditions.  We  ignore 
important  problems  such  as  fluid-structure  interaction,  non-rectilinear  motions,  very  low  Re  (below  10000), 
gust  and  other  unsteady  external  forcing  effects,  and  on  and  on. 

Having  established  something  of  a  baseline,  it  is  now  possible  to  be  more  systematic  about  pursuit  of  the 
various  problems  beyond  the  scope  of  the  present  study.  Before  passing  towards  enumeration  of  such 
possibilities,  we  pause  to  recapitulate  the  main  findings  from  Chapters  3-8. 


9.2  GENERAL  CONCLUSIONS 

At  the  juncture  between  fluid  mechanics  and  MAV  engineering,  we  are  interested  in  the  accurate  prediction 
of  aerodynamic  forces  -  lift,  thrust/drag,  pitching  moment  and  so  forth  (presently  we  are  limiting  focus  on 
longitudinal  quantities),  as  functions  of  motion  kinematics,  Reynolds  number  and  model  geometry. 
The  flowfield  structure  is  important  in  so  far  as  we  connect  flowfield  features  such  as  LEV  formation, 
growth,  and  shedding,  with  variations  in  the  aerodynamic  force  time  history,  such  as  dynamic  stall  spikes  in 
lift  coefficient. 

In  this  study  we  have  considered  primarily  6  cases:  “deep  stall”  and  “shallow  stall”  of  three  configurations: 
a  SD7003  airfoil,  a  nominally  2D  flat  plate  of  2  -  3%  thickness  and  round  edges,  and  an  AR  =  2  flat  plate  of 
the  same  sectional  shape.  These  were  studied  primarily  at  a  reduced  frequency  of  0.25  in  sinusoidal  periodic 
established  motion.  Broadly,  we  note  that  for  the  deep  stall  cases  of  all  configurations,  it  is  generally  easier 
to  capture  the  lift  coefficient  time  history  than  for  the  shallow-stall  case.  Likewise,  deep-stall  flowfields  vary 
less  from  shape  to  shape  than  do  the  shallow-stall  cases.  Lift  coefficient  history  is  easier  to  capture, 
in  general,  than  drag  coefficient  history,  and  pitching  moment  is  the  most  ornery  of  all.  The  importance  of 
leading  edge  vortex  resolution  for  force  coefficient  prediction  goes  in  the  same  sequence. 
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Classical  analytical  methods  -  inviscid  quasi-steady  lift  curve  slope  and  Theodorsen’s  formula  -  deserve 
mention  as  predictions  of  lift  coefficient  time  history  for  preliminary  estimates,  for  example  in  parameter 
studies  or  conceptual  design.  We  generally  find  that  both  shallow  and  deep  dynamic  stall  at  least  in  part 
effectively  produce  a  trailing  edge  Kutta  condition  where  statically  the  airfoil  would  be  stalled.  For  the 
SD7003  in  pure -plunge  (deep  stall),  the  simple-minded  2*pi*alpha  has  an  error  of  -0.4  CL  at  both  peak 
and  trough,  while  Theodorsen’s  formula  has  an  error  of  -0.2;  this  should  be  regarded  in  the  context  of  a 
peak-to-peak  lift  coefficient  variation  of  about  2.5.  The  error  bounds  are  somewhat  larger  for  the  2D  flat- 
plate,  but  absence  of  experimental  data  renders  conclusions  tenuous.  On  the  other  hand,  2*pi*alpha  is  a 
surprisingly  successful  approximation  for  the  AR  =  2  plate,  in  both  pitch-plunge  and  pure -plunge. 
So  caution  is  required  in  applying  Theodorsen’s  formula  and  quasi-steady  attached- flow  theory,  and  of 
course  one  should  appreciate  the  bounds  of  studying  only  a  few  cases.  Further,  this  predictive  capability  is 
generally  suited  only  to  lift  -  not  drag  and  especially  not  pitch.  Resolving  the  flowfield  correctly  is 
important  for  calculating  aerodynamic  loads  that  depend  on  the  pressure  distribution  explicitly,  and  not 
just  on  its  integration.  Even  so,  the  apparent  suitability  of  the  classical  theoretical  prediction  for  the 
AR  =  2  plate  is  surprising  and  merits  a  rigorous  parameter  study  of  aspect  ratio,  effective  angle  of  attack 
and  other  parameters. 

Not  surprisingly,  the  physics  of  boundary  layer  transition  -  facility  flow  quality  and  computational 
turbulence  modelling,  for  example  -  become  less  important  in  overall  flowfield  and  aerodynamic  load 
prediction  when  there  is  a  mechanism  forcing  transition,  such  as  round  leading  edges  of  thin  plates. 
Thus  the  scatter  in  our  results  for  plate  cases  is  much  less  than  for  airfoil  cases. 

Lacking  prior  experience,  one  might  assert  that  inviscid  methods,  and  even  RANS,  should  predict 
aerodynamic  force  more  accurately  for  shallow  stall  than  for  deep  stall,  because  presumably  complex 
viscous  and  transitional  effects  are  more  prevalent  for  separated  flows.  Our  experience  shows  that  for  the 
infinite-span  cases  examined  here,  this  is  precisely  not  the  case.  Scatter  amongst  the  computations  and  the 
experiments  for  shallow-stall  was  larger  than  for  deep  stall,  and  again  larger  for  the  airfoil  than  for  the 
plate. 

Only  one  finite  aspect  ratio  case  was  considered:  AR  =  2.  From  these  results,  in  deep  and  shallow  stall, 
we  find  generally  a  smaller  separation  and  a  tighter  LEV  in  the  AR  =  2  case  than  in  the  wall-to-wall  (that  is, 
2D)  cases,  but  no  discemable  difference  in  LEV  stability.  This  should  be  considered  in  terms  of  insect-wing 
LEV  stability  mechanisms  widely  found  in  the  literature,  but  again  is  only  one  data  point  and  requires  both  a 
deeper  examination  and  a  wider  parameter  study. 

At  high  reduced  frequency,  on  the  other  hand,  Theodorsen’s  formula  and  vortex  particle  methods  are 
excellent  predictors  of  lift  coefficient  time  history,  because  lift  is  traceable  mathematically  if  not  physically 
to  so-called  non-circulatory  or  inertial  effects,  which  the  inviscid  theory  captures.  Of  course,  this  is  limited  to 
periodic  flows  in  periodic  acceleration. 

Finally,  we  made  a  brief  examination  of  comparison  between  longitudinal  motions  and  a  rotating  motion 
or  “swinging”  wing,  where  the  latter  has  3D  effects  such  as  spanwise  pressure  gradients.  We  did  not  find  a 
large  difference  in  vortex  formation  or  lift  time  history  between  translation  and  rotation,  at  the  Reynolds 
numbers  and  motion  parameters  of  the  present  study.  This,  perhaps,  is  surprising  in  the  context  of  recent 
literature  -  and  again  merits  both  deeper  study  and  a  larger  parameter  space  investigation  -  of  which  we 
make  suggestions  below. 


9.3  RECOMMENDATIONS  FOR  FURTHER  WORK 

Broadly,  we  suggest  the  following  avenues  for  further  research. 
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9.3.1  Extensions  of  Cases  Studied  in  This  Work 

The  first  category  is  parameter  studies  based  on  the  already  achieved  results,  in  particular  to  extend  the 
findings  on  applicability  of  quasi-steady  and  Theodorsen’s  model,  to  explore  the  limits  of  effective  angle 
of  attack,  to  study  aspect  ratio  effects  and  to  find  clear  cases  where  leading  edge  vortices  materially  affect 
the  lift  coefficient  time  history.  Ideas  include: 

•  Sweep  of  the  pitch  parameter  X;  variation  of  X  such  that  effective  angle  of  attack  is  held  nearly 
constant,  by  varying  h  while  keeping  k  constant,  or  vice  versa,  or  both;  there  are  many  possibilities, 
all  nominally  equivalent  if  quasi-steady  reasoning  is  valid.  The  purpose  is  to  assess,  for  a  given 
value  of  maximum  effective  angle  of  attack  attained  in  the  motion,  indeed  whether  effective  angle  of 
attack  makes  sense. 

•  Parameter  study  to  ascertain  how  common  it  is  for  there  to  simultaneously  be  a  moment-stall  but 
no  lift-stall  for  sinusoidal  pitch-plunge,  plunge  or  pitch  problems. 

•  Extension  to  higher  Strouhal  numbers  to  study  deeper  dynamic  stall,  to  see  if  finally  one  finds  a  lift- 
stall  and  strong  dynamic  lift  due  to  the  leading  edge  vortex  formation,  shedding  and  convection. 

•  Further  exploration  of  blockage  effects  and  spanwise  variations  in  ground-test  facilities,  for 
example  by  testing  the  same  sectional  geometry  and  Re  but  with  different  physical  span  to  chord 
ratios  (all  nominally  wall-to-wall)  or  by  running  computations  in  3D  that  include  the  tunnel  walls 
and  the  model  support  system. 

•  Parameter  study  of  aspect  ratio  effects;  AR<2  and  >2.  Ideally  this  would  involve  volumetric 
velocimetry  where  one  obtains  not  only  three  components  but  also  a  volume  of  data  that  includes  the 
mid-span  region  and  the  tips. 

•  Extension  of  the  recent  studies  on  flapping-wing  propulsive  efficiency  [18],  to  see  for  example  how 
the  various  k-h-X  combinations  affect  thrust  production  and  power  consumption.  This  includes, 
of  course,  more  robust  direct-force  measurements  in  experiments,  and  more  thorough  benchmarking 
of  computations  in  predicting  drag/thrust. 

•  Extension  of  waving-wing  problems  to  the  high  angles  of  attack  encountered  in  insect 
aerodynamics,  where  LEV  production  and  retention  is  most  critical  to  flight  performance. 

•  Extension  of  waving-wing  problems  to  much  lower  Reynolds  number,  bringing  propeller-like 
experiments  [100]  with  those  studied  under  the  present  work. 

9.3.2  Conceptual  Departures  Towards  Other  Kinds  of  Motions 

Here  we  suggest  possibilities  in  keeping  with  abstracted  problems,  but  considering  questions  beyond  those 
raised  in  the  present  study.  Examples  include: 

•  Transient  longitudinal  motions  such  as  the  classical  problem  of  impulsive  start  of  translational 
motion,  pitch  ramps  or  plunge  ramps,  and  other  non-periodic  motions  (trigonometric,  linear  or 
otherwise)  where  there  are  unsteady  effects  not  only  from  the  attained  rate  of  motion,  as  would  be 
the  case  with  periodic  problems  too,  but  also  from  the  start-up. 

•  Oscillatory  problems  in  hover,  where  the  nominal  reduced  frequency  is  infinite  (because  there  is 
no  free-stream)  and  vortices  shed  from  each  motion  period  remain  nearby,  to  affect  the  flowfield 
and  aerodynamic  force  history  in  the  subsequent  period. 

•  Combination  of  longitudinal-plane  and  non-longitudinal-plane  motions,  where  the  pitch  varies 
sinusoidally  but  the  model  is  rotating  about  a  fixed  point  on  one  wing  tip. 
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9.3.3  Passage  to  More  Flight-Like  Problems 

Ultimately  the  purpose  of  MAV-inspired  aerodynamics  research  is  to  build  better  MAVs,  on  a  foundation  of 
scientific  knowledge  rather  than  mere  intuition.  To  achieve  this,  the  scientific  community  must  transition  its 
work  to  more  flight-like  problems,  and  to  abstract  questions  encountered  in  flight  applications  to  problems 
amenable  to  rigorous  study  in  the  lab.  This  has,  of  course,  been  a  problem  in  all  aeronautical  applications, 
but  the  very  novelty  of  MAVs  makes  the  old  problem  new  again.  Examples  of  more  flight-like  investigations 
include: 

•  Study  of  passive  fluid-structure  interaction  for  flexible  flight  surfaces.  This  is  a  huge  area,  but  one 
example  generalizing  our  work  would  be  flexible  plates  that  passively  camber,  akin  to  bird  wings, 
to  attenuate  leading  edge  separation  at  high  instantaneous  effective  angles  of  attack,  and  hopefully 
produce  better  propulsive  efficiency  from  flapping. 

•  Adaptation  of  both  periodic  and  non-periodic  motions  to  model  atmospheric  gusts  relevant  to  MAV 
flight  in  urban  environments.  This  involves  both  understanding  gust  spectra  (spatial  and  temporal), 
and  how  to  realize  them  in  a  ground  test  facility,  for  example  by  shuttering  the  tunnel  or  by  using 
the  motion  rig  to  impose  a  motion  (at  high  rate  these  two  are  demonstrably  not  equivalent). 

•  Study  of  motions  in  the  frequency,  amplitude  and  range  of  degrees  of  freedom  relevant  to  aggressive 
manoeuvres  of  natural  flyers,  such  as  insects  and  bats. 

Of  course,  no  such  list  can  ever  presume  to  be  comprehensive.  Our  only  ambition  is  for  it  to  be  cogent  and 
relevant  -  and,  sadly,  ineluctable  biases  stemming  from  one’s  own  work  make  such  an  assessment 
impossible. 

But  let  us  conclude  with  optimism:  Micro  Air  Vehicle  aerodynamics  is  fascinating  in  its  complexity, 
but  this  complexity  need  not  be  daunting  for  practical  engineering  approaches  to  aerodynamic  force 
modelling  for  vehicle  design.  We  already  have  a  basic  understanding  of  how  flapping-wing  MAVs  and 
even  natural  flyers  operate,  or  should  operate.  What  remains  is  a  thorough  assessment,  with  the  rigor  and 
enthusiasm  that  aeronautical  engineers  applied  to  fixed-wing  manned  aircraft  in  the  1920s  and  1930s, 
now  recast  in  the  new  century  of  powered  flight. 
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